o NANA NA NA 


TOU OI 


Revista de Engenharia 


Publicação da Associação dos Es: 
to Instituto Sunerior Técnico 


Técnica - revista de Engenharia 


Propriedade: Associação dos Estudantes do 
Instituto Superior Técnico 
Av. Rovisco Pais 1000 Lisboa Portugal 


Director Científico: Manuel de Abreu Faro 
Director Executivo: Paulo David Simões 
Coordenação Editorial: Elias Azevedo 


Coordenação Administrativa: Gerardo Vieira Lisboa 


Direcção, Administração e Publicidade: 
Associação dos Estudantes do 
Instituto Superior Técnico 
Av. Rovisco Pais 1000 Lisboa Portugal 
Telefones directos: 8481018 -8489323 
Extensão IST: 1248-1260 


Registo DGCS 101407 
Depósito Legal 41683 
ISSN 0040-1714 


A, mt 5000 exemplares 

Publicação Trimestral | 

Preço de capa:300$00 assinatura anual (Portugal): 1000$00 
(Estrangeiro): USD$12 


Contribuiram para este número: 
Junta Nacional de Investigação Científica e Tecnológica 
Conselho Directivo do Instituto Superior Técnico 
Impressão: ARTE 2, Lda. 


Técnica - número único de 1989 Impressão: Dezembro de 1990 
EDITORIAL | o 2 
Artigo Convidado 

SOBRE A EXCITAÇÃO DE ONDAS COMETÁRIAS HIDROMAGNÉTICAS 


Armando Larcher Brinca 3 
Notas Técnico-Científicas 


SOBRE A EQUIVALÊNCIA DAS FORMULAÇÕES DA ELECTRODINÂMICA DOS 
MEIOS EM MOVIMENTO 


Manuel de Abreu Faro 9 
POTÊNCIA E ATENUAÇÃO EM GUIAS DE ONDA 
João Luís Sobrinho 17 


EFFECT OF CURRENT ON GOLD ELECTROWINNING WITH FLUIDIZED BED 

ELECTRODE 

C.A.€. Sequeira and F.D.S. Marques 21 
MODELAÇÃO FÍSICO-MATEMÁTICA DOS FENÔMENOS OCORRENTES DENTRO 

DE CÂMARAS DE COMBUSTÃO 


Maria da Graça Carvalho 27 
SPECTRAL ESTIMATION IN THREE STEPS 
M.D. Ortigueira and J.M. Tribolet “43 


Técnica - 89 


EDITORIAL 


Num passado não distante a TÉCNICA configurou-se numa revista plena 
de interesse, integrando trabalhos de vanguarda nos domínios da ciência e 


da técnica. 
Espera-se que a partir de 1991 retome o ritmo perdido. 


Caminhar para esses objectivos e dar realidade às grandes finalidades de 
uma revista com a tradição da TÉCNICA é tarefa que apela 
necessariamente para os docentes, investigadores e alunos do 1.5.T. e bem 
assim para todos aqueles que cultivam os domínios científicos e técnicos 
que servem a Engenharia. 

Para relançar a revista há que acreditar na TÉCNICA, na acção que pode 
desenvolver na comunidade técnico-científica e, nomeadamente, ter 
presente que, em princípio, os artigos da TÉCNICA são referidos nos 


Abstracts das diferentes especialidades científicas. Assim tem acontecido. 
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RESUMO 


Uma breve introdução ao estudo dos cometas, e au 
seu papel na evolução do planeta Terra, precede a 
análise de instabilidades de ondas hidromagnéticas 
estimuladas por iões cometários de hidrogénio e 
oxigénio. Utilização de diagramas de Brillouin 
generalizados clarifica a natureza física dos modos 
ondulatórios e mostra que a interacção entre os 
dois tipos de iões é geralmente fraca, com cada feixe 
excitando independentemente (em quase todo o 
plano de dispersão) instabilidades ressonantes 
homólogas. O ambiente cometário sugere os 
parâmetros do modelo adoptado, mas os resultados 


são utilizáveis na interpretação de outras' 


observações espaciais de actividade ondulatória 
hidromagnética. 


Precedemos a descrição das ondas cometárias 
hidromagnéticas de uma breve introdução 
motivadora do estudo dos cometas. Estes corpos são 
antepassados do Homem no sentido lato, e poderão 
ter contribuído decisivamente para o seu 
aparecimento. Ao investigarmos os cometas, 
investigamos as nossas origens. 


PRIMÓRDIOS 


Para o sistema solar, tudo terá começado há 4,6 
milhares de milhões de anos numa nébula de gás e 
poeira em rotação. Por conservação da quantidade de 
momento angular, a sua contracção gravitacional é 
acompanhada de um aumento da velocidade de 
rotação. A interacção das forças de gravidade e 
centrífuga conduz à formação de um disco de 
acrescimento. 

A condensação central eventualmente ultrapassa o 
limiar da ignição nuclear e origina o Sol. 
Instabilidades gravitacionais na nébula solar 
precipitam a formação rápida de inúmeros objectos 
com dimensões lineares da ordem de alguns 
quilómetros.Estes protocometas, por colisões e 


Comunicação apresentada à Classe de Ciências da Academia 
das Ciências de Lisboa, na sessão de 21 de Abril de 1985. 


ABSTRACT 


A brief introduction to comets and their role in the 
evolution of our planet precedes the study of 
hydromagnetic wave instabilities fed by coexisting 
newborn (cometary) hydrogen and oxygen ions. 
Utilization of generalized Brillouin diagrams 
clarifies the physical nature of the modes and 
shows that the interaction between the two ions is 
generally weak, with each beam exciting resonant 
instabilities without undue influence from the 
other newborn ion species in most dispersion 
domains. Although cometary environments 
suggest the parameters of the adopted model, the 
results are helpful in the interpretation of other 
hydromagnatic wave activity in space. 


interacções gravitacionais, (il) evoluiram para os 
planetas, e (ii) foram expulsos (por influência de 
Jupiter e Saturno) do sistema solar, ou (iii) injectados 
(por influência de Urano e Neptuno) na Nuvem de 
Oort, zona exterior ao sistema solar que se estende 
até distâncias heliocêntricas da ordem de 1,6 anos 
luz. 

Os cometas foram assim criados dentro da nébula 
solar, um pouco antes da formação dos planetas, há 
4,6 milhares de milhões de anos. Tanto os cometas 
como os planetas são agregados de matéria 
interestelar condensada. A diferença reside na 
intensa reestruturação física e química a que foram 
submetidos os planetas durante a sua evolução, 
enquanto os cometas, conservados em condições 
ideais na Nuvem de Oort, se mantiveram 
relativamente imunes à acção do tempo. - 
Durante a translação do sistema solar em torno do 
centro galáctico da Via Láctea, perturbações 
gravitacionais estelares podem alterar algumas das 
órbitas (quasi-circulares) dos cometas na Nuvem de 
Oort para órbitas elípticas muito excêntricas que, por 
interacção com planetas, eventualmente se 
convertem em órbitas elípticas de períodos curtos 
(inferiores a 200 anos). Estas visitas cometárias 
transportam, portanto, uma fonte de matéria pristina 
e indicativa da composição da nébula solar. 
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Sobre a Excitação de Ondas Cometárias Hidromagnéticas 


Colisões cometárias 


As visitas dos cometas ao sistema solar interior 
podem terminar por colisão com o Sol ou planetas. No 
caso da Terra, reconhecendo que os cometas são 
basicamente bolas de neve sujas [Whipple, 1950] com 
traços significativos de várias substâncias como o 
metano e a amónia, as colisões têm sido invocadas em 
vários contextos. 

Parte importante da água nos oceanos poderá ter 
origem cometária. Também, se a atmosfera primitiva 
da Terra era pobre em hidrogénio (predominância de 
azoto e anidrido carbónico), os cometas podem ter 
transportado os ingredientes necessários à evolução 
prebiótica. Decomposição do vapor de água, metano e 
amónia por acção da radiação solar ultravioleta ou 
de descargas atmosféricas, seguida de recombinação, 
pode originar formaldeido e ácido cianídrico; estes, 
nas águas amoniacais da Terra primitiva, podem 
recombinar para formar os aminoácidos mais 
simples, ingredientes básicos das proteinas. 
Identicamente, poder-se-iam formar os componentes 
essenciais dos ácidos nucleicos. Estaria assim 
assegurada a síntese de moléculas orgânicas 
necessárias ao aparecimento da vida. 

Mas nem sempre as colisões com cometas estão 
associadas à origem da vida. Evidência crescente, 
mas ainda controversa, aponta para uma colisão 


importante há cerca de 65 milhões de anos, na - 


transição do período geológico Cretáceo para o 
Terciário. O inverno prolongado originado no 


bloqueamento da luz solar, com o consequente 


arrefecimento e interrupção da fotossíntese, poderá 
ter sido responsável pelas extinções maciças de 
espécies vegetais e animais nesta época, incluindo o 
desaparecimento dos dinoussauros. Curiosamente 
esta evolução abriu caminho à gradual 
predominância dos mamíferos e, espera-se, poderá 
contribuir para a sua sobrevivência por acção 
pedagógica dissuasora: poderiam ser idênticas as 
consequências de um inverno nuclear (provocado 
pelo Homem) e de um inverno cometário prolongado 
(cuja responsabilidade não cabe aos dinossauros). 
(Alguns investigadores julgam ter encontrado uma 
periodicidade da ordem dos 28 milhões de anos nas 
extinções maciças ocorridas na geo-história. 
Possíveis explicações para a intensificação das 
colisões cometárias com essa periodicidade incluem 
(1) a invocação de uma estrela companheira do Sol, a 
Nemesis, com uma órbita associada a esse período 
que, ao atravessar a Nuvem de Oort, origina as 
necessárias perturbações gravitacionais, e (ii) o 
movimento oscilatório do sistema solar 
relativamente ao plano da nossa galáxia, também 
com um período dessa ordem de grandeza.) 


Lendas, mitos e Halley 


Não surpreenderá que o aparecimento, por vezes 
exibindo formas espectaculares, dos cometas nos céus 


nocturnos dos nossos antepassados os tenham 
fortemente intrigado. Coincidência fortuitas das suas 
passagens com alguns acontecimentos importantes 
alimentaram uma mitologia, com tónica catastrófica, 
que só gradualmente sucumbiu ao progresso dos 


conhecimentos. 


A referência mais antiga a cometas chegada aos 
nossos dias é uma frase chinesa de 1500 anos a.C.: 
“Quando (o imperador) Chieh executou os seus fiéis 
conselheiros, apareceu um cometa". A tapeçaria de 
Bayeux sugere que o cometa de 1066 foi responsável 
pela deposição do rei Harold de Inglaterra, imposta 
pelos invasores Normandos. Montezuma II do México 
considera que a chegada de Cortez foi anunciada por 
um ominoso cometa. À observação por Halley, na sua 
infância (8-9 anos), dos cometas de 1664, associado à 
Grande Peste de Londres, e de 1665, coincidente com 
o Grande Incêndio de Londres, poderá ter marcado os 
seus interesses científicos. 

Edmond Halley, nascido em 1656, publicou o seu 
primeiro artigo aos 18 anos no "Philosophical 
Transactions" da "Royal Society of London"; 
intitulava-se "Um método directo e geométrico para a 
determinação dos afélios, excentricidades e 
proporções (das órbitas) dos planetas primários que 
dispensa a hipótese da igualdade do movimento 
angular". Com 21 anos dirigiu uma missão à ilha de 
Santa Helena para efectuar o levantamento das 
estrelas do hemisfério celeste austral. Em 1684 
desloca-se a Cambridge para discutir com Isaac 
Newton a órbita dos planetas. Newton envia-lhe 
porteriormente o manuscrito em latim "De motu 
corporum in gyrum” onde prova que a lei da atracção 
universal com o inverso do quadrado da distância 
implica a satisfação das três leis de Kepler. 
Apercebendo-se da importância deste trabalho, 
Halley regressa a Cambridge e insiste com Newton 
para expandir as suas ideias num livro. Halley acaba 
por pagar a publicação dos "Principia Mathematica" 
porque a "Royal Society” não dispunha, na altura, de 
verbas para o efeito. 

Em 1705 Halley edita "Uma sinopse da astronomia 
dos cometas" onde, por aplicação das Leis de Newton, 
determina as órbitas de cometas e postula o 
reaparecimento do cometa que hoje tem o seu nome 
(curiosamente, Newton inclinava-se para trajectórias 
cometárias parabólicas). De sublinhar, no entanto, 
que a influência das contribuições científicas de 
Halley resistiria à amputação deste seu trabalho. De 
facto, foi ainda um fundador da geofísica, com 
especial interesse pelo geomagnetismo e o seu 
ecletismo pode-se inferir da publicação em 1691, aos 
35 anos, de artigos sobre: a distância Terra-Sol, o 
comportamento do vapor de água na atmosfera, a 
matemática das grandezas infinitas, a espessura de 
filmes de ouro, a invasão de Inglaterra por Júlio 
César, e a História Natural de Plínio. Escreveu o seu 
último artigo científico, sobre um eclipse lunar, aos 
81 anos, quatro anos antes de falecer. 
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Cometologia moderna 


A desmitificação resultante das investigações de 
Halley introduziu o estudo dos cometas na sua era 
moderna. Ao interesse pelo cálculo das órbitas, 
sucederá até ao presente a curiosidade pelo cometa 
propriamente dito (sua origem e constituição), pela 
sua interacção com o meio ambiente, e pelo seu papel 
na evolução terrestre. 

A seguir à Segunda Guerra Mundial, Ludwig 
Biermann conclui, por observação das caudas dos 
cometas, que deveria existir um fluxo permanente de 
partículas originado no Sol, o que veio a ser 
confirmado directamente pela nave soviética Luna 3, 
em 1959, com a primeira observação do vento solar. 
Há cerca de dez anos, um grupo de cientistas da 
Universidade da Califórnia em Berkeley que incluia 
Luis Alvarez, prémio Nobel de Física e seu filho 
géologo Walter Alvarez, investigou a constituição de 
uma fina camada de barro cinzento separando 
sedimentos correspondentes aos períodos Cretáceo e 
Terciário, em Gubbio, Itália. Verificou-se que esse 
barro era significativamente (cerca de 30 vezes) mais 
rico em irídio do que as camadas contíguas. Esta 
constatação, posteriormente confirmada noutros 
pontos do globo, sugere a colisão da Terra com um 
corpo extraterrestre há cerca de 65 milhões de anos, 
aquando da extinção maciça de espécies animais e 
vegetais, na transição do Cretáceo para o Terciário. 
Isto porque os espectros estelares e a constituição dos 
meteoritos mostram que a matéria extraterrestre é 
mais rica em vários elementos, entre os quais o 
irídio. (Tendo as matérias terrestre e extraterrestre 
uma origem interestelar comum, as suas diferentes 
composições parecem paradoxiais, mas são 
justificáveis. Os elementos em que a matéria 
extraterrestre é mais rica -- irídio, ósmio, ródio, ouro, 
platina, -- tendem a acompanhar, quando fundidos, o 
ferro; na recém formada Terra, com as rochas ainda 
fluidas, o ferro concentrou-se no núcleo líquido 
central, arrastando com ele estes elementos e 
criando, portanto, a correspondente rarefacção 
superficial.) Mais recentemente, variações 
detectadas na razão das concentrações dos isótopos 
87 e B6 de estrôncio para a transição sedimentar em 
causa sugerem a ocorrência de intensa chuva ácida 
que poderá ter contribuído para a extinção de 
espécies. 

Com a aproximação do regresso do cometa Halley em 
1986, já em plena época de exploração espacial, 
vários grupos de investigação aproveitaram a 
oportunidade para realizar pela primeira vez 
observações “in situ”. Os europeus, em homenagem 
ao autor do fresco 'Adoração dos Magos' (Capela de 
Arena, em Pádua), onde o cometa de Halley na sua 
aparição de 1301 é utilizado como estrela de Belem, 
designaram a missão espacial por Giotto (di 
Bondone). Os japoneses enviaram duas sondas, 
Sakigake ("pioneiro") e Suisei ("cometa"). Os 
soviéticos lançaram as naves Vega 1 e 2 


Armando Larcher Brinca 


("VEnera + GAlley" = Venus + Halley; a língua russa 
não tem a letra H, usando o G aspirado em sua 
substituição) que passaram primeiro por Venus. A 
missão ISEE-3 foi adaptada à missão ICE 
("International Cometary Explorer"), por forma a 
observar o cometa Giacobini-Zinner e, a grande 
distância, o cometa Halley. No ano anterior, tendo 
em vista uma melhor compreensão da fenomenologia 
associada à interacção do vento solar com cometas, a 
missão AMPTE ("Active Magnetospheric Particle 
Tracer Explorers”) tinha criado 'cometas artificiais' 
através da libertação de iões de bário na bainha 
magnética (região do geoplasma entre o choque em 
arco e a magnetopausa). 

De todas estas missões resultou uma impressionante 
acumulação de dados (de que as ondas cometárias 
hidromagnéticas representam uma componente 
muito restrita), objecto de aturada análise e 
interpretação, e, não surpreendentemente, um 
intensificado interesse científico na realização de 
missões cometárias mais ambiciosas. Está neste caso 
o projecto CRAF ("Comet Rendezvous and Asteroid 
Flyby"), a lançar em Agosto de 1995; trata-se da 
primeira missão do tipo Mariner Mark II, sendo 
curioso referir estar conjuntamente programada a 
missão Cassini*, a Saturno e ao seu satélite Titan 
cuja termodinâmica foi objecto de comunicação (4 de 
Fevereiro de 1988) à Academia das Ciências de 
Lisboa pelo Professor Jorge Calado. 

A missão CRAF, para além de observar o asteróide 
Hamburga, em contraste com as missões pontuais ao 
cometa Halley, faz o acompanhamento do cometa 
Kopff desde antes do afélio até depois do periélio. 
Haverá assim aportunidade de analisar em detalhe a 
evolução do cometa na sua aproximação do Sol e 
concomitante interacção com o vento solar, e 
determinar a sua composição com a utilização de 
uma sonda que penetrará o núcleo do cometa. 


Ondas cometárias: observações e ressonância 
ciclotrónica 


As observações de turbulência hidromagnética 
efectuadas pelas missões cometárias são dignas de 
realce. Por exemplo [Tsurutani e Smith, 1986], a 
turbulência que envolve Giacobini-Zinner atinge 
distâncias cometocêntricas da ordem do milhão de 
quilómetros, e intensidades três ordens de grandeza 
acima da actividade ondulatória média no vento 
solar. No que se segue, tentaremos identificar fontes 
de energia livre e instabilidades ondulatórias 


* O astrónomo italiano Gian Domenico Cassini foi nomeado por 
Luis XIV como primeiro director do Observatório de Paris. 
Descobriu quatro satélites de Saturno e o primeiro indício 
(Divisão de Cassini) da existência de estrutura nos anéis deste 
planeta, Sugeriu a Halley que o grande cometa de 16B0 era o 
mesmo que tinha sido observado por Tychu Brahe em 1577; 
embora errada, a sugestão implicava a periodicidade do 
movimento de (alguns) cometas e pode ter desempenhado um 
papel seminal no trabalho cometário de Halley. 
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potencialmente responsáveis pela excitação desta 
turbulência hidromagnética ("Hidromagnético" 
significa aqui terem estas ondas electromagnéticas 
frequências menores do que a frequência híbrida 
inferior do hidrogénio, no referencial do vento solar). 
Ao aproximar-se do Sol, a superfície do cometa 
sublima-se. Os gases libertados são eventualmente 
ionizados pela radiação ultravioleta, ou por troca de 
carga eléctrica com partículas do vento solar. Os iões 
cometários recém criados apresentam uma 
velocidade instantânea relativamente ao vento solar 
de —400 km/s, e a trajectória que traçam depende da 
orientação relativa (ângulo a) da velocidade do vento 
solar, V., € do campo magnético interplanetário, 
B,. No referencial do vento solar, a componente 
paralela (a B,) da velocidade iónica, V., cos a, 
permanece constante, enquanto que à componente 
perpendicular, V., sina, está associada a um 
movimento ciclotrónico circular. Forma-se no espaço 
das velocidades uma distribuição em anel 
perpendicular, animada de um movimento paralelo 
de deriva. Surgem assim duas fontes de energia 
iónica que, em condições adequadas, poderão 
estimular a excitação de modos ondulatórios: a 
energia perpendicular e a energia paralela de deriva. 
Um mecanismo que torna possivel a transferência de 
energia das partículas cometárias para as ondas, é a 
ressonância ciclotrónica. Num magnetoplasma, as 
ondas electromagnéticas paralelas têm 
necessáriamente polarização circular (esquerda ou 
direita, conforme rodam no sentido do movimento 
ciclotrónico das cargas positivas ou negativas, 
respectivamente), rodando os seus campos num ponto 
do espaço com a velocidade angular correspondente à 
sua frequência w. Cargas eléctricas sem velocidade 
paralela exibem um movimento circular puro com 
velocidade angular 9 diferente, em geral, de w. 
Ressonância ciclotrónica (e, portanto, forte 
interacção onda-partícula) ocorrerá se, (i) existindo 
movimento paralelo das particulas com velocidade V, 
o associado efeito de Doppler kV (onde k é o número 
de onda da excitação com frequência w) fizer 
coincidir a frequência ciclotrónica da partícula com a 
frequência da onda por ela sentida, OQ = w-kV e, (ii) 
os sentidos de rotação dos campos da onda e da 
velocidade perpendicular da partícula forem 
coincidentes. 

A visualização das condições de ressonância é 
facilitada pelo traçado de diagramas de Brillouin 
onde se representa a dispersão das ondas e o efeito de 
Doppler no feixe de partículas. A representação 
clássica num quadrante único, com w e k positivos, 
torna-se, no entanto, confusa por requerer 
consideração simultânea da polarização das ondas, 
sinais das cargas eléctricas, e sentido relativo de 
propagação. Mais atraente [Brinca e Tsurutani, 
1988] é a utilização dos quatro quadrantes do plano 
de Brillouin, onde ocorre naturalmente a separação 
dos modos em função da sua polarização e sentido de 
propagação em relação à velocidade paralela de 


deriva das particulas ressonantes. Todas as 
intersecções ocorridas no plano de Brillouin entre as 
curvas de dispersão dos modos e do feixe 
correspondem a ressonâncias ciclotrónicas reais. 
Como decorre da topologia das ondas 
hidromagnéticas paralelas [Brinca e Tsurutani, 
1988], uma grande variedade de interacções 
onda-partíicula pode ter lugar, muito embora seja 
impossível (por motivos fisicamente óbvios) a 
ressonância ciclotrônica entre feixes positivos e 
ondas de polarização circular direita com sentidos de 
propagação opostos. 

Frisamos que, apesar destes diagramas ajudarem à 
interpretação física da interacção, a globabilidade da 
ressonância ciclotrónica não fica caracterizada 
completamente. Em particular, (i) não se infere 
quando a interacção origina crescimento, ou 
atenuação, da onda, (il) a energia perpendicular do 
feixe de partículas não é considerada, e (iii) existem 
instabilidades não ressonantes cuja natureza 
dispensa a ocorrência da interacção ciclotrónica. Na 
particularização para os cometas devemos ainda ter 
presente que (a) existem vários tipos de iões 
cometários (por exemplo, 10es de hidrogénio e 
oxigénio); (b) a velocidade de deriva paralela dos iões 
cometários é dada por V.. cosa; (c) as grandes gamas 
de frequência e número de onda a cobrir para o estudo 
das ondas hidromagnéticas aconselha a adopção de 
escalas logarítmicas nos diagramas de Brillouin, o 
que vai distorcer as formas familiares das curvas de 
dispersão dos feixes e das ondas. 


Ondas cometárias: Modelação e análise linear 


O modelo adoptado no estudo da estabilidade dos 
modos hidromagnéticos num ambiente cometário 
utiliza o vento solar como referencial. Assume-se 
composto por populações maxwellianas de protões e 
electrões, com densidade comum N,= IN = 95/em” e 
temperaturas T =2x10º K e T =8x10º K; a sua 
velocidade relativamente ao o 1 tem um módulo 
V..=400 km/s e o campo magnético interplanetário, 
com uma amplitude B =8 nT, forma um ângulo a 
(parâmetro variável do modelo) com V |. Este vento 
solar é permeado por populações cometárias ténues 
de iões de hidrogénio (p)), oxigénio (o) e 
correspondentes fotoelectrões (e), com densidades 
N=N=N,/2=0,05/cm”. As distribuições de 
velstidado al são separáveis em V e V| (orientações 
relativas a B), a distribuição paralela é 
maxwelliana com uma temperatura de 0,01 eV e 
velocidade de deriva V cosa, e a distribuição 
perpendicular, de forma arbitrária (não afecta a 
estabilidade dos modos de propagação paralela), 
apresenta um valor quadrático médio V * sinº a. O 
estado de equilibrio deste magnetoplasma com cinco 
populações tem, portanto, carga global e corrente 
nulas. 

A linearização das equações de Maxwell e Vlasov 
para amplitudes complexas dos campos variando 
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como exp i (K.r-wt) define a equação cinética de 
dispersão D(w=w +iy, K) = 0 do.meio em análise. 
Num problema de valores iniciais, o vector de onda K 
é real e alinhado com o eixo dos zz (e com B ); modos 
instáveis têm y>0. As partículas cometárias 
apresentam uma velocidade de deriva ao longo de B , 
V .. cos a, e os modos electromagnéticos paralelos com 
polarização circular (esquerda ou direita) 
propagam-se no mesmo sentido, ou sentido oposto, da 
deriva cometária. O programa WHAMP [Rônmark, 
1982] resolve numericamente a equação de dispersão 
exacta, e os resultados são apresentados por forma a 
caracterizar a estabilidade hidromagnética do meio 
em função do ângulo a, e das temperaturas e 
densidades cometárias [Brinca e Tsurutani, 1988]. 


Ondas cometárias: Resultados e conclusões 


A Figura 1 mostra os resultados obtidos para a= 45º, 
valor típico associado à espiral do campo magnético 
interplanetário a uma distância heliocêntrica de 
uma unidade astronómica. O painel superior mostra 
o diagrama de Brillouin; as taxas de crescimento dos 
modos predominantemente associados com os feixes 
de protões e oxigénio ionizado estão indicadas nos 
paineis médio e inferior, respectivamente. Os modos 
são identificados pela sua polarização circular (L - 
esquerda, R - direita), sentido de propagação (F - com 
a deriva dos feixes cometários, B - contra), feixe 
cometário (P - protões, O - oxigénio) e frequência (H - 
“alta”, S - baixa”). 

Cada feixe excita instabilidades idênticas. Os modos 
LFPH e LFOH são ondas ciclotrónicas próprias dos 
feixes: a sua ocorrência não depênde da existência do 
vento solar de fundo. Estas oscilações (ondas 
ionico-ciclotrónicas) com w=Q (s=p, 0) no 
referencial dos feixes, satisfazem, por efeito de 
Doppler, a w=S2 + kV cosa no referencial do vento 
solar. À sua fonte de energia livre reside na 
anisotropia (excesso de energia perpendicular 
relativamente à energia térmica paralela) dos iões 
cometários. Os modos RFPH e RFOH exibem taxas 
de crescimento máximo na vizinhança das 
interacções das curvas de dispersão dos modos RF do 
vento solar não perturbade e das curvas de dispersão 
dos feixes, confirmando a sua natureza ressonante 
ciclotrónica. Os feixes cometários impõem ainda o 
acoplamento entre os modos LB (LBPS e LBOS) e os 
modos RF (RFPS e RFOS), criando os modos mistos 
MP e MO. A transição entre os modos LB e RF 
implica a ocorrência de estruturas não oscilatórias 
(«o =0), puramente crescentes (y >0). 

Como discutido em pormenor por Brinca e Tsurutani 
[1988], a topologia das instabilidades 
hidromagnéticas é alterada quando ocorrem 
modificações no vento solar de fundo (variações do 
ângulo a, em particular), ou nas populações 
cometárias (densidades e temperaturas). Análise da 
evolução das características dos modos 
intervenientes permite extrair conclusões relevantes 
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Figura 1 

(Painel superior) Diagrama de Brillouin generalizado dos modos 
electromagnéticos paralelos excitados no modelo nominal do 
ambiente cometário, para Q = m/4. (Paineis inferiores) Taxas de 
crescimento correspondentes. Às curvas a tracejado representam 
as curvas de dispersão dos feixes. A notação utilizada na 
identificação das instabilidades está definida no texto. Os 
quadrantes superiores, w,>0, (inferiores, w,<0) do plano de 
Brillouin estão associados a modos com polarização circular 
esquerda (direita). Nos quadrantes kw,>0 ( kw,<0) os modos 
propagam-se em sentido idêntico (oposto) ao da deriva paralela 
dos iões cometários. A frequência angular ciclotrónica e o raio de 
giração térmico dos protões do vento solar normalizam as escalas 
axiais, O, = 211 x 0,122 rad's e pp = 474km. 


para a interpretação das observações da actividade 
ondulatória em ambientes cometários: 

- É razoável estimar os efeitos dos dois feixes 
cometários (iões de hidrogénio e oxigénio) como 
resultando da sobreposição da dispersão de cada um 
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(a sua interacção mútua é fraca). De notar, no 
entanto, que a coexistência de feixes cometários de 
iões pesados com massas similares origina 
fenomenologia distinta [Brinca e Tsurutani, 1989 b). 

- Os modos LFPH e LFOH não suportam aumentos 
(mesmo moderados) na temperatura paralelos dos 
feixes. 

- Os feixes de protões geram taxas de crescimento 
ondulatório mais fortes do que os feixes de oxigénio 
quando as densidades cometárias são elevadas. À 
medida que nos afastamos do núcleo do cometa (e 
estas densidades decrescem), a hierarquia do 
crescimento altera-se e os modos associados aos iões 
de oxigénio predominam (marginalmente). 

- À alteração introduzida na dispersão básica do 
magnetoplasma com a adição de iões pesados 
(oxigénio) origina o aparecimento de um novo modo 
electromagnético com polarização circular esquerda 
cujas propriedades de propagação oblíqua (polariza- 
ção, taxa de crescimento e compressibilidade) 
merecem, pela sua especificidade, investigação 
separada [Brinca e Tsurutani, 1987]. 

A aplicação dos resultados desta análise sobre o 
comportamento das instabilidades hidromagnéticas 
paralelas à turbulência cometária e outros 
ambientes espaciais deve ainda ser complementada 
com investigações similares para propagação oblíqua 
[Brinca e Tsurutani, 1989 a], e com o estudo 
(possivelmente com códigos de simulação nimérica) 
das correspondentes evoluções não lineares dos 
modos instáveis. 
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RESUMO 


Apresentam-se quatro das mais conhecidas 
formulações da Electrodinâmica dos Meios em 
Movimento: a de Minkowski, a de Chu, a 
Amperiana e a de Boffi. 

Como Tai tem vindo a demonstrar, são 
equivalentes. 

Na nossa opinião a de Minkowski prima pela 
simplicidade, e por isso a adoptámos. 

Seguimos um método de demonstração que 
complementa a demonstração de Tai, e com isso 
julgamos ter contribuido para reforçar a convicção 
da equivalência das diversas formulações. 


1. INTRODUÇÃO 


Foi sob o título de "Sobre a Electrodinâmica dos 
Corpos em Movimento" que Einstein, em 1905, 
apresentou a Teoria da Relatividade Restrita [1]. 
No entanto, não entrou no pormenor das equações 
referentes a meios materiais [2). 
Minkowski, em 1908, em plena posse do “princípio de 
relatividade" foi quem finalmente resolveu o 
problema da Electrodinâmica dos Meios em 
Movimento. 
Assim, ou quase assim, nos diz sobre o assunto 
Sommerfeld [2]. 
Fundamentalmente, trata-se de uma teoria 
macroscópica estruturada a partir de quatro 
entidades E, B, De H, duas a duas constituindo 
tensores anti-simétricos de segunda ordem: 
Tensor campo electromagnético -—[E,B|>F. 

(1) 


Tensor excitação =>[D, H] > Gº 


* Actividade desenvolvida no âmbito do Centro de 
Electrodinâmica, da UTL, do INIC. 

Comunicação apresentada à Classe de Ciências em 17 de Julho 
de 1986. Memórias da Academia das Ciências de Lisboa, Tomo 
XXVII, 1986. 


ABSTRACT 


Four of the most widely known formulations of 
Electrodynamics in Moving Media are presented: 
Minkowski's, Chu's, Boffi's and Amperian 
formulation. 

In accordance with Tai's recent demonstration, 
they are equivalent. 

In our opinion, Minkowski's formulation stands out 
for its simplicity, and we have, therefore, adopted 
it. 

Having followed a demonstration method that 
complements Tai's demonstration, we believe to 
have contributed to strengthen the assumption of 
the equivalence of the various formulations. 


É notável que as equações de Maxwell 


aB 
VXE =-— V.B=0 
dt 


— àD 
VxH=J+ — V.D=p (2) 
dé 
consintam essa formulação tensorial. 
As equações de Maxwell são compatíveis com a 
equação de continuidade associada ao movimento das 
cargas eléctricas 


dp 
Vd+ — =0 (3) 
dt 


Não se tratando de uma teoria arquitectada a partir 
de distribuições físicas de cargas e correntes era 
natural que se viesse a defrontar com outras. 

Assim sucedeu. 

Das diversas formulações possíveis, [3], apenas 
trataremos de quatro o que é suficiente para elucidar 
a natureza da questão abordada. 

A formulação de Minkowski não recorre a modelos. 
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Às outras teorias recorrem a modelos, mas isso não 
implica que fisicamente o modelo exista. Basta que 
permita calcular correctamente o campo 
macroscópico. 

Como Tai tem vindo a demonstrar, [3] e [4], as 
diferentes formulações são equivalentes. 

O mesmo concluem Penfield e Haus, |5). 

Como tantas vezes tem sucedido, apenas haverá 
razões de simplicidade suficiente, de elegância 
conseguida, em última análise apenas haverá lugar 
para uma opção essencialmente estética. 

Apesar de ser assim, não estamos certos de que a 
questão se tenha encerrado. 

Usando uma aproximação de natureza física 
demonstraremos a referida equivalência. 


2. QUATRO POSSÍVEIS FORMULAÇÕES 
Até hoje salientaram-se as seguintes formulações: 


M - Minkowski, 1908 [Sommerfeld, 1952]. 
Entidades fundamentais [E,B] e [D,H!]: 


Teoria EBDH. 
Entidades definidas 
P=D- E, E (4a) 
B 
— — E] (4h) 
H 


Não recorre a modelos. 


Chu, 1960. 

Entidades fundamentais [E,H] e [P,M): 
Teoria EHPMv. 

Recorre a modelos: 

Uma distribuição de dipolos eléctricos 


P=nql 
(5) 
Uma distribuição de dipolos magnéticos 


u,M E Mio 


Explicita a velocidade v do meio em relação ao 
referencial de repouso macroscópico. 

Como é evidente, n en são densidades volúmicas, q, 
e q, cargas eléctricas e magnéticas e | el as 
distâncias orientadas entre cargas. 

A - Formulação Amperiana. 

Panofsky e Philips, M = 0, 1950, 

Fano, Chu e Adler, 1960. 

Entidades fundamentais [E,B] e [P,M]: 
Teoria EBPMv. 


Recorre a: 
Uma distribuição de dipolos eléctricos 


p= nal. 
(6) 
Uma distribuição de correntes eléctricas 


M=nla 
a E 


Em que n representa uma densidade volúmica de 
correntes eléctricas l. circulando em espiras 
elementares a que correspondem áreas orientadas a . 
Explicita a velocidade relativa ao repouso, v. 


B - Boff, 1957. 
Entidades fundamentais [E,B] e [P,M]: 
Teoria EBPM. 
Na base das distribuições da formulação 
Amperiana: 


P* = nql 
| (7) 
M* = nla, 


utiliza uma polarização e uma magnetização 
equivalentes definidas por 


(8) 
Mé=Mº+PÊxu 


3. COMPARAÇÃO PRÉVIA DAS QUATRO 
FORMULAÇÕES 


No referencial onde a matéria está em repouso os 
campos de idêntica designação coincidem, 
observando-se as seguintes correspondências: 


M- Minkowski E B D-c,E (B/y,)-H 
C- Chu E mo(H + M) P M 
(Ja) 
A - Amperiana E B g M 
B - Boffi E B P M 
Verificam-se, portanto, as relações 
P= D-z,E M=(B/g,)-H (9b) 


No vácuo as quatro teorias reduzem-se a uma teoria 
comum que está bem comprovada. 

Em qualquer das teorias, J e p constituem um 
quadrivector 


(Ji)=[J,cp] (9c) 


10 Técnica - 89 


Os tensores momento-energia conduzem aos mesmos 
efeitos mecânicos embora se ofereçam com diferente 
partição. 


4. FONTES E VÓRTICES DO CAMPO E 
CORRESPONDENTES EQUAÇÕES DE 
MAXWELL 


Na base dos modelos descritos na secção anterior, 
podemos escrever imediatamente as fontes e vórtices 
das equações de Maxwell. 

Basta aplicar as equações da Electrodinâmica no 
vácuo a cada uma das distribuições consideradas e 
ter presente as definições das densidades de 
correntes eléctrica e magnética de deslocamento. 
Deste modo se obtêm correntes e cargas cujas 
densidades se unem formalmente aos outros vórtices 
e outras fontes de campo. 

Na formulação de Minkowski não há outras 
distribuições que não sejam pe J e assim teremos: 


Fontes do campo: 
p>D 0—B (10a) 
Vórtices do campo: 
ab) | 9B 
(J,—)5H (—)5E (10b) 
dt dt 


As equações de Maxwell têm uma escrita covariante 
que é: 


aD 
(M) didi dad: 7.D= p 
(10c) 
VxE=-— V.B=0 
Na formulação de Chu: 
São fontes do campo: 
(0,-V.uM>n H (p,-VP)- E (1la) 


São vórtices do campo: 


oa UM 
tu—,——.. VxuMxW>E 
“mx dt R 
(11b) 
dP 
Jde—,;— VxPxwy>H 
O dt di 


Em conformidade, as equações de Maxwell 
escrevem-se: 
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e V.E=pV.P 


dE dP 
VxH=J+e— + — + VxPXxuv 
o dt dt 
(11c) 


(C) 


au, H au M 


VxE=- -— .V Xu MXu 
di ÚU 


pVH=ViuM) 


Admitiremos que esta escrita é covariante qualquer 
que seja a velocidade relativa ao repouso, v. 


Na formulação Amperiana: 


São fontes do campo: 
MXvuv 
(p,-V.P,Y. E J>e E 0—B (12a) 
c 
São vórtices do campo: 
aB 
(- ria 
(12b) 
dk aP à MXv. B 
de — — VxPXv, VXM,-— | — 
od dl dt H, 


Nestes termos, as equações de Maxwell escrevem-se: 


dk aP | d MXu 
-—VYxB=J+te—+— +VxPxvu+VxM-— —— 
Ho om à à 
Mxvuv 
e VE=p-VP+YV. 
á c 
(A) (12c) 
aB 
VYXE=-— 
ol 
V.B=0 
Admitiremos que esta escrita é covariante. 
Na formulação de Boff: 
São fontes do campo: 
(p,-VP)>e É 0 —>B (13a) 
São vórtices do campo: 
aB 
-— |->E (136) 
dt 
dk qP B 
de —,— VXxM)>» — 
“dt di 
[é 
1d, 
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Vejam-se as relações (8) entre (P,M) na formulação 
Amperiana e os (P,M) utilizados por Boffi. 
As equações de Maxwell escrevem-se: 


l dk  àP 
Ny XD=) de pos + TAM 
o di dl 


u 
(B) lá (13c) 
e VE =p-V.P 
aB 
VXE=-— 
dl 
VB=0 


Mais uma vez se admitirá a covariância da escrita 
das equações de Maxwell. 


5. INTERPRETAÇÃO FÍSICA DAS FONTES E 
VÓRTICES DOS CAMPOS 


Na distribuição de Chu, e como é bem conhecido, [6], 


p,=-VP (14a) 
traduz uma densidade volúmica correspondente a 
uma polarização não uniforme no espaço. 

A esta polarização está ainda associada uma 
densidade de corrente 


aP 
Jp = (14h) 
ol 
e uma corrente de convecção 
J =ppu=-vV.P=VxPxu (14c) 


Pe 


O mesmo se dirá para a distribuição de dipolos 
magnéticos: 


=-V.M Jdua= 


oM 
Mo Ma 


J v=-vV.M=VxMXu (15) 


Mc PM 
Na formulação Amperiana teremos, como 
anteriormente, 


aP | | 
Pp =-V.P Jdo=— Jo =VxPxXvuv 


P di Pe 


Da magnetização resulta um termo devido à não 
uniformidade de M no espaço, |6], 


dy = VxM (16h) 
e ainda 
Mxu 7 M x sap (iéo 
p = v = = à ADC 
MR e? e? 

o que significa uma polarização 
= Mx uv ed 
Puro 2 dd 


que surge e apenas se pode explicar por um efeito 
relativista, [6]. 

Se atendermos a que a corrente que circula na espira 
é parte de um quadrivector, imediatamente se 
conclui que, em movimento, uma espira, neutra no 
referencial de repouso, aparece com uma distribuição 
polarizada de carga eléctrica. 

Claro está que esta polarização, além de contribuir 
para a densidade volúmica (16c), contribui também 
para a densidade de corrente 


J Bão (16) 
=. — ? je 
PR al e? 


Na formulação de Boffi duas considerações são feitas: 
A densidade de corrente de convecção 


py=VxXPxv 


define uma magnetização equivalente 


M,=PXv (17a) 
que podemos adicionar a M. 
De modo idêntico a polarização 
po MXu 
ui = 2 (17) 


se poderá adicionar a P. 
Deste modo a polarização P* e a magnetização M* da 
formulação Amperiana conduzem às relações (8) 


PE pe. Mé=Mº+PAxu 


Mº x uv 
2 


Cc 


(16a) o que faz desaparecer v das equações de Maxwell. 
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6. TRANSFORMAÇÕES DAS ENTIDADES 
ELECTROMAGNÉTICAS NAS DIFERENTES 
FORMULAÇÕES 


Consideremos três referenciais de inércia, S,5S'e S”. 
Admitiremos que a matéria está em repouso 
macroscópico no referencial S. 

S está animado de velocidade v em relação a S'. 

S” está animado de velocidade u em relação a S'. 

A transformação de Lorentz, 

o facto de [J, cp] constituir um quadrivector; 

a covariância, postulada, das equações de Maxwell, 
permitem-nos obter as leis de transformação dos 
campos, de S' para S" e reciprocamente. 

Admitindo assim, obtém-se: 


Formulação de Minkowski: 
E"=E' +y(E, +uxB') 
DB" + vB -u X Et /e') 


(M) 


D"'=D' +y(D! +u x M'ye) 
E" = M', + vai -«uxD') 


Formulação de Boffi: 
E" 

Idênticas às da formulação de Minkowski 
B" 


(B) 
P"=P' + y(P-ux MV e) 
Mº = M', + Yi + uxP') 
Formulação de Chu: 
(18) 
Ee.” + y(E, +uxyg H) 
Hº =H + Y(H! -uxe E) 
(C) 


=" our, +yu xXx (vx Pe 
Mº=M +yM +yux(vxM)/CM" 


Formulação Amperiana: 
E"=E +yYE +uxBº) 


B"=B' +y(B)-ux E//c*) 
(A) 
P” 
Idênticas às da formulação de Chu 
M" 
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7. OBTENÇÃO DAS RELAÇÕES ENTRE AS 
ENTIDADES INTERVENIENTES NAS 
DIVERSAS FORMULAÇÕES 


No estado de repouso as entidades com a mesma 
designação coincidem ou relacionam-se por: 


D=c E+P B=u(H+M) (19) 


No caso anterior podemos fazer coincidir S com S" de 
onde resulta u=v. 

S" será então o referencial onde a matéria está em 
repouso macroscópico. 

Os campos em S” são assim idênticos em todas as 
formulações ou correspondem-se através de (19). 

Nas relações (18) faremos, então, u = v e em seguida 
substitui-se v por (-v) a fim de representar F' em 
função de F”. 

Os campos F" serão todos indicados pelo mesmo 
símbolo. 

Os campos F' designam-se por F*, Fº, Fº ou Fº 
consoante se trate da formulação de Minkowski, de 
Chu, Amperiana ou de Boff. 

Assim procedendo obtém-se: 


E“=E +y(E -v x B) 

Bº=B +y(B + vX E //c) 

D“=D,+y(D,-v x H//cº) 

Hº=H +y(H +vxD) 

E*=E +y(E -v x B) 

Bº=B +y(B + vxX E /€) 

PP=P+y(P +vxM) 

MB=M +y(M -v XP) 

(20) 

EC=E +ytE -v xp H) 

HC=H +ylH +vxe E) 

PC=P +yP' 

MC=M +yM, 

E*=E, + v(E -v x B) 

Bº=B, + y(B + vx E//€º) 
PP=P +yP 

MA=M, + yM 
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De (20) facilmente se obtém o seguinte quadro de 
equivalências: 


EDBH EPBM 
Minkowski Boffi 
EM E* 
DM *g EP + Pp? 
BM -B 
HM Bº/u - Mº 
EPHMyv EPBM 
Chu Amperiana 
ES-v x uM E“ 
ER +P coB*+P* + vx Mº/e 
p (He + M9 Bº 
H'+v x Pº B*/u -M*-v x P* 


Como Tai salienta, [3], [7], no trabalho de Fano, Chu 
e Adler (1960); (Ec + wo v x He) e (Hc - e, v x Ec) são 
introduzidas e postuladas como sendo as forças 
associadas às cargas unitárias, eléctrica e magnética. 
Mas do quadro de equivalências resulta facilmente: 
EM+v x B“=E +y vx H' 

(21) 
HH y x Dº=H". Ev E“ 


e então compreende-se bem que a força resultante é a 
mesma em qualquer das teorias, a de Minkowski e a 
de Chu. 

Este um aspecto interessante e muitos outros se 
poderiam considerar. À conclusão é sempre a mesma: 
as formulações são equivalentes. 


8. SOBRE O METODO UTILIZADO 


O método utilizado consistiu em admitir que as 
equações de Maxwell dizem sempre o mesmo na 
relação que estabelecem entre o campo e as fontes e 
os vórtices desse campo que será [E,B] ou [D,H). 

As fontes contribuem para a divergência dos campos 
e os vórtices para os rotacionais. 

Nas formulações que não as de Minkowski, as fontes 
e vórtices unem-se formalmente, em todos os efeitos 
físicos previstos, a pe J e, bem assim, às densidades 
de correntes de deslocamento. 

Então, os campos calculam-se utilizando a 
electrodinâmica no vácuo: é como se todo o espaço do 
meio material se convertesse no vácuo de parâmetros 
constitutivos [u, e ] e nesse vácuo existissem, de 
facto, as distribuições de carga e de corrente (as 
verdadeiras e as equivalentes). 

Foi esta atitude que nos permitiu escrever 
imediatamente na Secção-4 as equações de Maxwell. 
Para que haja consistência é necessário que, na base 
do quadro de equivalências estabelecido, as equações 
de Maxwell escritas numa dada formulação resultem 


das equações de Maxwell escritas nas outras 
formulações. 

Consideremos então as equações de Maxwell na 
formulação de Minkowski: 


an” 
vxHM=J+ E VD” =p 
(22) 
as” 
vx E! = a vB” =0 


Se utilizarmos, por exemplo, a formulação de Chu 
obtém-se das equivalências estabelecidas: 


d 
VX(EvXuM)=-—p (UH +MO) Vie E +PO=p 
| (23) 
| Í 
Vx(H +uxP )=J+ ag (eg + PO) Vu (H+M)=0 
Mas estas são exactamente as equações obtidas por 
Fano, Chu e Adler (1960) e facilmente se convertem 
nas equações (11lc) escritas na Secção 4 


relativamente à formulação de Chu. 
O mesmo se passa com todas as outras formulações. 


9. CONCLUSÕES 


Demonstrámos a equivalência formal das 
formulações da electrodinâmica dos meios materiais 
em movimento. 

Diferença física, para se evidenciar, exige 
experiência em meios materiais em movimento onde 
é difícil ou até impossível medir os campos. 

Esta a argumentação que normalmente se levanta. 
Mas, como se verificou, nas "forças de Lorentz" 
intervêm os campos [E”, B”] ou [E*, H'], (formulação 
de Minkowski e formulação de Chu). 

Desde que a medição de um campo envolva um efeito 
mecânico é natural que seja dificil dizer se existe uma 
formulação mais correcta do que as outras. 

- | Notrabalho apresentado, obtivemos exactamente 
os resultados já anunciados por Tai [3, 4] e por 
Penfield e Haus [5]. 

No entanto, o método que se seguiu é diferente. Por 
exemplo, Tai, consegue a equivalência das equações 
de Maxwell do que se segue o quadro de 
equivalências. 

Nós, escrevemos imediatamente as equações e, em 
seguida, utilizando as transformações que garantem 
a covariância das equações de Maxwell e o facto de 
em repouso os campos coincidirem nas diversas 
formulações, obtivemos o quadro de equivalência. 
Depois, a partir desse quadro, mostramos que as 
equações de Maxwell, nas diversas formulações, 
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podem deduzir-se umas das outras. Nisto consiste a 
originalidade deste trabalho. 

No que respeita a Penfield e Haus as metodologias 
também são essencialmente diversas. Trata-se da 
obra que mais informação contém sobre a 
Electrodinâmica dos Meios em Movimento. 
Salienta-se isto para eventual consulta. 

- Note-se, isso é importante, que os efeitos 
relativistas evidenciados na Secção 5, relações (16c) 
e (16e) conjugados com a equivalência de todas as 
formulações evidenciam bem que as equações de 
Maxwell são essencialmente relativistas. 

- Finalmente pode perguntar-se: qual a melhor 
formulação? 

Em face do exposto, é talvez uma pergunta sem 
sentido. 

A de Minkowski foi a primeira e é formalmente a 
mais simples e elegante e isso afigura-se-nos 
importante. 
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POTÊNCIA E ATENUAÇÃO EM GUIAS DE ONDA 
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SUMÁRIO 

Neste trabalho, o autor analisa de uma forma 
simples e original, julgamos, a expressão da 
potência num guia de ondas, associada a modos TM 
ou TE, quando expressa na componente 
longitudinal do campo. A argumentação utilizada 
permite ir além, esclarecendo sobre o 
comportamento assimptótico da constante de 
atenuação, dos referidos modos, quando a 
frequência tende para infinito. 


1- INTRODUÇÃO 

A determinação da potência (P) e atenuação (a), 
aproximação de primeira ordem, num guia de ondas 
envolve, para cada modo de propagação, o cálculo de 
integrais da forma: 


| 
- À | HP. HT ds 
2 P C 


2P 


As expressões obtidas são exaustivas mas não 
evidenciam o comportamento em frequência de P ea. 


No presente trabalho mostra-se que, a partir de 
propriedades gerais do campo electromagnético num 
guia de ondas, se pode de imediato prever a evolução 
de Pe a com a frequência. As conclusões a que se 
chega são independentes da geometria da secção 
transversal do guia. 

A figura 1 esclarece o significado dos simbolos 
intervenientes no texto. 


2- CAMPOS DOS MODOS TM E TE 


Muito sumariamente pode dizer-se que as 
propriedades gerais dos campos, num guia de ondas, 
se consubstanciam nas relações que abaixo se 
escrevem para uso ulterior (consultar qualquer livro 


ABSTRACT 


In this work, the author analyses in a simple and 
original way the formula representing the power in 
a waveguide, associated with TM or TE modes, 
when expressed on the axial component of the field. 
The line of reasoning used makes it possible to 
predict the asymptotic behaviour of the 
attenuation constant, as the frequency approaches 
infinity. 


figura 1 


sobre guias de onda, por exemplo Propagação Guiada 
do Prof. Abreu Faro). 
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Modos TM 


la) A resolução do problema de Dirichlet 


V2E +kºE =0 1a) 
L E CU É 


[8,=0) 2a) 


conduz a ÉE,. 


lla) As componentes transversais dos campos são 
dadas por: 


kKSVE 
ET= (É) tz 3a) 
k k 
E =] E) id 
k. k Z La 
Modos TE 
Ib) A resolução do problema de Von Neumann 
Vv2H +KkH =0 1) 
L ZE CL 
(= a o) 2b) 
dn 


C 


conduza H,. 


Ilb) As componentes transversais dos campos são 
dadas por: 


HT = -5( + | 3b) 
ko! k 
L t 
: ko VH Xe : ab) 
E - (E) " & =lmptl xe, 
C 


em que km e Zm São, respectivamente, a constante de 
propagação e a impedância característica de onda do 
meio no interior do guia, que se admite sem perdas 
(om = 0). | 


1b), 3b) e 4b) são equações duais de la),3a) e 4a) o que 
justifica por estas derivarem, e apenas, das equações 
de Maxwell na ausência de fontes. Note-se que não 
estão envolvidas as equações 2a) e 2b) que traduzem 
condições fronteira em superfícies condutoras 


perfeitas. 


3- POTÊNCIA 


É bem conhecida a fórmula da potência, associada a 
modos TM ou TE, expressa na componente 
longitudinal do campo: 


Modos TM 
2 
1 1 í a 1 | —* 
p=>2>(>) vi (=) | E E'as, >? 
2 LO f gs 22 T 
m c T 
Modos TE 


ps t. + ; 
p= ZhT) 1(<) R H Hds, 


em que f. é a frequência crítica do modo de 
propagação que se considere. 

Como se mostrará estas expressões são fáceis de 
apreender. Evidenciar este facto é o primeiro 
objectivo do trabalho. 

Verifica-se de Ila) e IIb) que, para cada solução do 
problema de Dirichlet, modos TM, ou do problema de 
Von Neumann, modos TE, as componentes 
transversais do campo crescem com a frequência. A 
componente transversal de nome contrário ao da 
componente longitudinal, isto é HT nos modos TM e 
ET nos modos TE, varia com 


e)=fg) 


À componente transversal do mesmo nome da 
componente longitudinal, isto é E! nos modos TM e 
HT nos modos TE varia com 


CE) e) =(EM (E) 


dr 
f 


A potência, dada de um modo geral por 


1 | 

Re | | (EX H* e dSy)=- | (ET x HT. e ds 

2 ]s. zo PF 2]s - e 
T T 

é proporcional ao produto das componentes 

transversais do campo. Vimos como cada uma destas 


evoluía com a frequência. Podemos então escrever 


pe(LJ A (E) 6) 
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A relação anterior foi estabelecida para cada solução 
de la) ou Ib). Não é portanto estranho que na fórmula 
da potência apareça, a multiplicar, um integral 
envolvendo a componente longitudinal do campo. 
Este integral depende da geometria da secção 
transversal do guia. Finalmente o factor 1/Zm, modos 
TM, ou Zm, modos TE, é necessário para que as 
expressões 5 a) e 5 b) fiquem dimensionalmente 
correctas. 


4- ATENUAÇÃO 


A cada modo de propagação está associado uma 
constante de atenuação, a , que caracteriza as perdas 
de potência nas paredes condutoras do guia 
(aproximação de primeira ordem). Como se sabe a é 
dado por: 


a=— 7) 
2P 
em que: 
1 = do 
P=-R H'. H' ds 8) 
P 2PJç 
=+| (Ex HO) e ds, 9) 


Rp é a parte real da impedância característica de 
onda das paredes metálicas. Recorda-se, desde já, que 
Rp «MV ff. 

Na secção anterior estudou-se a dependência do 
denominador de 7) com a frequência. A mesma 
atitude para com o numerador permite prever o 
comportamento assimptótico da constante de 
atenuação quando f tende para infinito. 


Modos TM 


Nestes modos, por definição, HZP = Q ao longo do 
contorno €C. Por outro lado HTP = 0 em C. (O autor 
não conhece uma demonstração da universalidade 
desta condição). 

A potência de perdas na parede é dada por: 


1 
P = =R,| HE. WE do 10 a) 
C 


P 2 


Atendendo a 4 a) é fácil verificar que: 


ee) MME) 1º 


C Lh C 


João Luis Sobrinho 


Conjugando 6)e 11 a) conclui-se que, nos modos TM, 
a constante de atenuação torna-se proporcional a Vf/ 
f. quando f tende para infinito. 


Modos TE 


Nestes modos HZP = () no contorno C. 
Se HTP = 0, ao longo de C, então a potência de perdas 
nas paredes é dada por: 


1 | 
Pp => | Hº”. Hds 10 b) 
p PJc 


2 
/£ 
Pa vy— 
pf 


Constata-se que: 
E 


11b) 


Conjugando 6) e 11 b) conclui-se que a constante de 
atenuação tende para zero quando f tende para 
infinito. Um exemplo deste caso, são os modos TE, 
nos guias de secção circular. 

Se, por outro lado, HTP = () em, pelo menos, um troço 
de C, então a 10 b) temos de adicionar o termo: 


E) EM=(5) | ME) 


o que torna a constante de atenuação proporcional a 
Vf/ f. quando f tende para infinito. 
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ABSTRACT 


Electrowinning of gold from caustic cyanide 
solutions made to simulate carbon column strip 
solutions was studied using a laboratory size 
side-by-side fluidized bed electrode (FBE). The 
variable considered was the current through the cell. 
“The results show that there are three competing 
cathodic reactions: gold deposition, oxygen reduction 
and hydrogen production. Oxygen reduction occurs 


because dissolved oxygen comes from the anode ' 


where it is produced. It is most significant at low 
currents. Hydrogen production occurs when the 
concentration of gold is low enough or the current is 
high enough that the gold diffusion limiting current 
is attained or nearly attained. 


INTRODUCTION 


Interest in gold production has risen substantially in 
the last few years. This is due to its relatively high 
price on the market while those of other metals 
remain low. Many medium and small size gold 
mining and milling facilities have opened in the past 
ten years. Most of the new plants use carbon 
adsorption techniques for purification and 
concentration of leach solutions. After the carbon is 
loaded with gold, it is stripped with a hot caustic 
cyanide solution. The gold in the resultant solution 
can either be electrowon or precipitated by 
cementation using zinc dust. The latter process is 
known as the Merrill-Crowe process. Because there 
is no need for deaeration or clarification of the 
solution with electrowinning, as with the 
Merrill-Crowe process, electrowinning is now a 
common procedure for recovering gold from strip 
solutions. Also there is interest recently in 
electrowinning of very dilute leach solutions. 

The same basic types of electrolytic cells have been 
used for the past 30 years to commercially electrowin 
gold from carbon strip solutions i. e., steel wool 


packed bed electodes. New electrowinning cells have 
been developed which can effectively electrowin gold 
from dilute solutions. The fluidized bed electrode 
(FBE) is one of these cells and it has been extensively 
tested since its invention in 1966. It has been tried 
for electrowinning of copper, nickel, zinc and other 
metals. These studies show that the fluidized bed 
electrode performs well in the electrowinning of 
metals from dilute solutions. 

One reason for the good performance of the fluidized 
bed electrode is that with a particulate cathode there 
is a high surface area which means metals can be 
won from dilute solutions at a low actual current 
density in a relatively small cell. Another advantage 
is the high turbulence in the bed due to fluidization. 
This and particle collision lead to a lowering of the 
diffusion layer thickness which correspondingly 
increases the limiting current for electrowinning. 
The goal of this study is to look at the reactions and 
other parameters involved in the electrowinning of 
gold from cyanide solutions with a fluidized bed 
electrode. In particular, the effect of current through 
the FBE cell is reported in this paper. 


EXPERIMENTAL 


The fluidized bed electrode used in this study is the 
side by side type of cell arrangement and is the same 
type used in most recent investigations of the FBE [1- 
4). The cell was made from welded polypropylene and 
had a fluidized cathode of conducting particles. These 
particles were relatively spherical copper shot, 20x30 
mesh, U.S. series. They were electroplated with gold 
and, therefore, behaved electrochemically like gold. 
The cathode had a volume of 50 cm3 when fluidized 
and contained 280 gm of gold coated copper particles 
which had an approximate surface area of 2700 cm2. 
The same particles were used for all experiments 
with the FBE since the diameter of the particles 
changed only about 2% total for all experiments 
combined. The bed was fluidized by flow of 
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electrolyte up through the cathode. A 100 mesh, U.S. 
series, stainless steel screen was used as the cathode 
feeder electrode and was placed opposite the 
diaphragm in the cathode compartment. It was about 
the same size as the diaphragm. 

The anode and cathode compartments were 
separated from each other with a polypropylene cloth 
diaphragm. A small portion of the electrolyte flowed 
from the cathode through the diaphragm to the 
anode compartment to form the anolyte. This 
electrolyte returned to the electrolyte from the 
cathode compartment at the top of the cell, and this 
mixed stream then entered a 2 É reservoir. The anode 
was made of stainless steel mesh and was flush 
against the diaphragm. There was no oxidation of 
the anode. 

In addition to the FBE cell, a two liter reservoir, a 
hot water bath, a power supply, a corrosion resistant 
pump, and two multimeters, were used in the 
experimental set-up. 

All solutions were made from doubly distilled water 
and reagent grade KAu (CN). 2H50, NaOH (1,0%) 
and NaCN (0.1%). The mixed electrolyte was placed 
in the two liter holding flask which was then put into 
the hot water bath. Once the solution reached the 
desired temperature it was sparged with nitrogen for 
two minutes to remove excess oxygen. Flow of 
electrolyte to the cell was initiated and adjusted to 
obtain proper fluidization (FBE) or flowrate. 
Samples of the electrolyte were taken at the same 
time that power to the cell was provided. Two 
samples were taken every 1.5 minutes thereafter, 
one from the holding flask and one from the top of the 
cell. The samples were analyzed for gold content 
using an Instruments Laboratory 351 Atomic 
Absorption Spectrophotometer. Each experiment 
was run until the solution was depleted to less than 1 
ppm of gold.Power was applied galvanostatically and 
the cell voltage was recorded with every sample 
taken. 

Concentration analysis of the samples from the 
holding flask as a function of time was used to 
construct a plot of gold concentration vs. time, Figure 
1, and the slope of this plot was used to calculate the 
current efficiency (CE) for gold deposition using the 
equation: 


CE(%)= 100(FV/A ud AuJ/dt 
=0.8166(V/Dd[ Au] (dt (1) 


where d/Au]/ dt is the slope of the gold concentration 
vs. time plot in ppm/min, V is the volume of the 
solution in liters, 1 is the current, F is the Faraday 
constant and AAy is the atomic weight of gold. 


Power consumption (PC) was also calculated using 
the cell voltage (E) with the equation: 


Power consumption, KWh/kg Au 
0 03 0.6 09 12 Ea 


Cell voltage, V 
[o 16 


m Concentration 

O Current efficiency 
150h 6 Cell voltage 

Y Power consumption 


IZ0 


90h 


Time, min 


60 


30 


100 ; 
Gold concentration, ppm 
' E ma E Si DST 
O 20 so 60 BO 0] 
Current efliciency, 


Fig. 1 Plot of raw data and calculated values as a function of time. 
PC(kWh/kg Au) = 8.333ElIMd[Au]/dt) (2) 


Both current efficiency and power consumption were 
calculated for every time interval that a sample was 
taken. 

The parameter varied in this study was the cell 
current, 0.2A to 2.0A i.e., 100 to 1,000 A/m2 
superficial current density for the FBE. The solution 
concentration and temperature were chosen to 
simulate commercial carbon column strip solutions. 


RESULTS AND DISCUSSION 


Raw data from each experiment was plotted in the 
format typified in Figure 1. The gold concentration 
vs. time curve is plotted from the analysis of samples 
taken from the 2 liter reservoir during the 
experiment. Cell voltage was measured to obtain the 
voltage vs. time curve. Calculated values of current 
efficiency and power consumption, using equations 1 
and 2, are also shown. The current efficiency initially 
increases due to equilibration of the entire system. 
Once equilibrium is reached the current efficiency 
starts to decrease as time proceeds. This is due to 
depletion of gold which also increases cell voltage. 


Gold Electrowinning Reactions 


In the caustic auro-cyanide system there are three 
basic reactions which occur at the cathode and 
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compete for the current. These reactions have been 
noted by others [5,6]. The most important and 
essential cathodic reaction is the plating of gold: 


Au(CNJ9 + é = Au + 2CN (3) 


The other two reactions, which reduce current 
efficiency, are evolution of hydrogen and reduction of 
oxygen: 


H$50 + e = 1/2H, + OH (4) 
1/20, + H50 + 2e = 20H (5) 


The predominant anodic reaction is the production of 
oxygen. This reaction provides the oxygen for 
reaction 5: 


20H" = 1/205 + H50 + 2e (6) 


This oxygen does not permeate the membrane into 
the cathode but rather is circulated through the 2 
liter reservoir to the cathode because of electrolyte 
recycle. Another anodic reaction of small magnitude 
is oxidation of the cyanide ion: 


20H" + CN = CNO + H50 + 2e (7) 


A minor reaction that takes place in the bulk 
solution is oxidation of the cyanide ion by dissolved 
oxigen: 


1/20, + CN' = CNO (8) 


Since reactions 3, 4 and 5 constitute all cathodic 
reactions, the current eficiency, CE, is given by 
CE(%) = 1001A4y/ (Lau + IH + Lo) (9) 
where lau, In, and Ig are the cathodic currents of 
reactions 3, 4 and 5, respectively. 
At the conditions of these tests (pH= 12.5 and 8000), 
the approximate Nernst potentials, Er, for the three 
respective cathodic half reactions, 3-5, with respect 
to the standard hydrogen electrode are as follows: 


Au(CN) 9/Au:Er=-0.60+ 0.070 log [Au(CN)9] / 


[CN-2]= 0.60V at 100 ppm (10) 
Ho0/Hs: Er = - 0.874V 

(11) 
02/0H": Er = 0.355 + 0.0174 log Po, (12) 


=0.33V 


The concentration dependent term of equation 12 
was evaluated by measuring the concentration of 
dissolved oxygen entering the cell, 1.6 ppm, and 
comparing it with the saturation level in equilibrium 
with oxygen gas, 21 ppm. 
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Current Efficiency 


The effect of current efficiency and power 
consumption is shown in Table 1. Power consumption 
will be discussed later. There are two features that 
stand out in the current efficiency curves and they 
are: 1) a pronounced maximum at 0.8A, i.e., 400 A/m2 
for the FBE,; and 2) increase in current efficiency 
with increased gold concentration. Since the current 
efficiency is given by Equation 9, it is necessary to 
consider all three cathodic reactions involved, Lau, 
ly, and Ig in terms of their respective currents to 
interpret data presented in Table 1. 

To better understand this, consider first two 
polarization curves for the cathode. The potential of 
the cathode was measured with respect to a standard 
calomel electrode by placing a Luggin capillary filled 
with electrolyte at the outer edge of the cathode near 
the diaphragm. These measurements, shown in 
Figure 2 are for solution with no gold (more negative 
curves) and solutions which contained greater than 
100 ppm of gold. Reactions 4 and 5, hydrogen 
evolution and oxygen reduction, are responsible for 
the shape of the curves with no gold. At low currents 
oxygen reduction consumes almost all of the current. 
When the current is at 0.25A i.e., approximately 
-1.24V (SCE), oxygen reduction lo, is at or near its 
limiting current, which can be seen by the change in 
the shape of these curves.The concentration of 
dissolved oxygen was measured to be 1.3 to 2.0 ppm. 
If all of the oxygen were consumed at the standard 
flow rate, 1.35 1/min, it would amount to a current of 
0.35 to 0.5A. Therefore, a limiting oxygen current of 
0.25A is of proper magnitude. Above about 0.25A 
(-1.24V) hydrogen starts to evolve, Ig, and consumes 
almost all additional current. 

With gold present the shape of the previous curve of 
Figure 2 changes. 

When the cathode has a potential above (less 
negative than) -1.1V (SCE) reduction of gold and 
oxygen occurs. At potentials more negative than 
about -1.1V, all three reactions take place 
simultaneously and compete for the current. An 
interesting point to note with the polarization 
curvesis that because of the variation of potential 
within three dimensional cathodes there are no well 
defined Tafel slopes and there is mixing of the 
different regions of the curves which causes the 
regional definition to be smoothed out. 

From the results of the polarization curves it is 
possible to explain the variation of the current 
efficiency vs. current, Table 1. At low current, left of 
the peak, the cathode is in the region of oxygen 
reduction. and gold deposition. Since oxygen 
reduction has a limiting current of 0.25 to 0.3A and 
also has a slightly positive Er as compared to the 
negative Er for the Au(CN)9/Au couple (see 
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solutions with no gold. 


equations 9 and 11), the 05/0H" couple will attain its 
limiting current before gold.deposition starts to take 
place i.e., low current efficiency for gold. But gold 
deposition increases with increasingly negative 
potential, while oxygen reduction stays at its same 
rate (its limiting current). Therefore, as total 
current increases there will be an increase in current 
efficiency as lay increases relative to Io. 

The polarization curves also show that hydrogen 
starts to evolve at just below the peak current of 0.8A 
or 400A/m? for the FBE (Table 1), while on the other 
hand, gold deposition reaches near its limiting 
current. Ás the potential becomes more negative 
hydrogen evolution, Iy, increases while both Ia, and 
lo stay the same which causes the decrease in the 
current efficiency to the right of the peak. 

The effect of Io and ly can be seen better if the data 
in Table 1 are plotted as gold current, Iau, vs. total 
current,l, (gold current = total current x CE). This is 
done in Figure 3 that shows that gold reaches its 
limiting current at low gold concentration. 

At higher gold concentrations gold deposition does 
not quite reach its limiting current i.e., lay is still 
increasing with increasing I. At very low currents no 
Ha is produced as previously explained and oxygen 
consumes 100% of the current. Figure 4 clarifies this 
and shows how the cathodic current is broken up into 
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Fig. 3 - Gold current as a function of total current at different 
gold concentrations. 


the three currents lay, ly and Lg for the FBE at high 
gold concentration. 

Another feature of the data on current efficiency vs. 
current, Table 1, is the increase in current efficiency 
with increased gold concentration. This occurs 
because at low gold concentrations there is mass 
transfer control of gold plating and mass transfer 
control depends on the gold concentration. 
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The power consumption data of Table 1 reflect mostly 
the trend of current efficiency (compare equations 1 
and 2) since the voltage increased relatively little 
with increasing current i. e., from about 1.5V to 2.5V. 
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CONCLUSIONS 


The following conclusions have been drawn from the 
data: 
1. Hydrogen evolution in the electrowinning of 
gold occurs when the concentration of gold is low 
enough or the current is high enough that the 
limiting gold diffusion current, lay, is attained or 
nearly attained. 


2. Dissolved oxygen consumes almost its full 
theoretical possible current (limiting current) in 
all experiments because of the high overpotential 
for oxygen reduction (equation 11) which is 
always attained in operation. This effect is seen 
most clearly in experiments with low current. 


3. A high electrolyte conductivity is needed to 
maintain a relatively high current efficiency in 
the FBE. 


It is clear that further evaluation of the cell is 
required, being necessary to study the effect of 
parameters like NaOH and NaCN concentrations, 
temperature of operation, flow rate, bed expansion, 
etc., so that optimised working conditions can be 
established for the electrowinning of gold with a 


fluidised bed electrode. These aspects will be treated 


in a future publication. 


Superficial 
Current Density,A/m2 | 
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Table 1 - Current efficiency and power consumption as a function of superficial 
current density 
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RESUMO 


Este artigo descreve o trabalho de modelação física 
e matemática e simulação numérica dos fenómenos 
ocorrentes dentro de câmaras de combustão 
desenvolvido no Departamento de Engenharia 
Mecânica do Instituto Superior Técnico. Os 
modelos físicos de turbulência, combustão, 
radiação e emissão de poluentes são descritos. O 
modelo matemático é constituído por um sistema 
de equações, diferenciais, não lineares às derivadas 
parciais, cuja solução passa pela discretização 
dessas equações, conduzindo a equações às 
diferenças finitas. O algoritmo foi aplicado ao 
cálculo de diferentes situações reais sendo 
apresentados os resultados numa demonstração 
das potencialidades do algoritmo. 


1. INTRODUÇÃO 
1.1 O Problema em Estudo 


No passado, o conhecimento empírico acumulado ao 
longo de dezenas de anos, juntamente com modelos 
teóricos aplicados na forma de balanços globais, 
determinaram a forma como era feito o projecto de 
câmaras de combustão e a procura das suas condições 
de funcionamento optimizadas. 

Apesar do desconhecimento da distribuição espacial 
das propriedades relevantes no interior das câmaras 
de combustão (i.e. velocidades, temperaturas, fluxos 
de calor e concentrações das espécies químicas), o seu 
funcionamento revelou-se economicamente eficiente 
até ao início da crise petrolífera. Contudo, a 
consciência da escassez das fontes de energia, o 
aumento do preço dos combustíveis e a legislação 
restritiva à emissão de poluentes libertados na 
combustão levaram ao aprofundamento do estudo dos 
fenómenos ocorrentes nos sistemas de combustão 
atrás referidos. 


ABSTRACT 


The present paper describes the work on 
mathematical modelling and computer simulation 
of the flow, heat transfer and combustion precesses 
in three-dimensional combustion chambers 
developed at the Mechanical Engineering 
Department of Instituto Superior Técnico. The 
physical models of turbulence, combustion, 
radiation and pollutant emissions are described. 
The mathematical model involves the solution of a 
set of partial differential. equations which are 
discretized using a finite difference technique. The 
algorithm was applied to several industrial 
situations and has proved to be an important tool 
for the furnace and boiler engineer designer. 


A simulação experimental em modelo real de novas 
câmaras de combustão é impensável devido aos seus 
proibitivos custos, ou mesmo, o estudo empírico de 
modificações na geometria de câmaras existentes, 
restando por conseguinte, a única alternativa viável: 
a modelação numérica de câmaras de combustão 
reais. À simulação em computador, através de um 
algoritmo de cálculo de uma câmara real que hoje é 
possivel, resultou da conjunção de esforços na 
modelação de fenómenos, tais como escoamentos 
turbulentos com ou sem reacção química e processos 
de transferência de calor. A resolução destes modelos 
matemáticos foi possível à custa dos avanços na 
análise numérica. O desenvolvimento nestas áreas, 
conjuntamente com o advento dos grandes 
computadores digitais, tornaram possível a solução 
de grande parte dos fenómenos práticos no domínio 
de Mecânica dos Fluidos, da Transferência de Calor e 
da Combustão. É hoje possível obter detalhes por 
exemplo, do escoamento e transferência de calor 
dentro de caldeiras , de câmaras de combustão de 
turbinas de gás ou de fornos industriais. 
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1.2 Revisão Bibliográfica 


A simulação numérica de câmaras de combustão 
tri-dimensionais foi iniciada por Patankar and 
Spalding [1,2]. Nestes trabalhos foi apresentada a 
base dum modelo numérico para escoamentos 
tri-dimensionais com combustão cujos principais 
pressupostos eram o uso de um novo modelo de 
combustão designado por "Sistema de Reacção 
Química Simplificada” (SRQS) e o cálculo da 
viscosidade efectiva através de uma relação 
algébrica. 

Na presente comunicação, a revisão bibliográfica vai 
ser dividida em três partes, baseando-se esta divisão 
no tipo de aplicação. 


1i) Câmaras de Combustão de Turbinas de Gás 


Patankar e Spalding [3] introduziram no modelo 
apresentado em [1,2] o então recente modelo de 
turbulência a duas equações desenvolvido por Jones 
and Launder [4]. Neste trabalho foram utilizadas 
coordenadas cartesianas e a câmara de combustão 
anular considerada foi idealizada como uma caixa de 
forma rectangular. 

Jones et al. [5] and Jones and Priddin [6] 
desenvolveram um modelo para a simulação dos 
fenómenos ocorrentes dentro de uma câmara de 
combustão de turbina de gás. O modelo permite 
tratar combustíveis líquidos, descrevendo o campo de 
concentração das potas através de uma equação para 
a função densidade de probabilidade do tamanho das 
gotas, e calcular a concentração de NO, assumindo 
que a formação de NO pode ser descrita pelo 
mecanismo de Zeldovich. O modelo de combustão 
empregue baseou-se no modelo de equilibrio químico 
de Gordon e McBride [7]. 


Serag-Eldin e Spalding |8] calcularam o escoamento, 
com combustão, numa câmara de combustão 
tri-dimensional, usando coordenadas polares e 
condições de fronteira cíclicas para o tratamento do 
"swirl”". O modelo de combustão utilizado foi o SRQS. 


Jones e McGuirk [9] usaram uma modelação análoga 


à de [6] para testar os resukados obtidos, mediante 
confronto com dados experimentais, para geometrias 
bi-dimensionais axissimétricas e tri-dimensionais, 
com ou sem combustão. Para câmaras 
bi-dimensionais axissimétricas, com combustível 
gasoso o acordo entre os resultados numéricos e 
experimentais é razoável, embora um conhecimento 
imperfeito das condições de fronteira torne difícil 
conclusões precisas e haja indicações de incapacidade 


do modelo de turbulência utilizado, modelo a duas ' 


equações, de prever com suficiente precisão o 
escoamento, particularmente quando há "swirl”. 


Boysan et al. [10] apresentaram um modelo para a 
combustão de sprays em câmaras de combustão de 
turbinas de gás. O modelo de combustão utilizado 
envolve a resolução de equações de conservação das 
fracções mássicas de combustível e de oxigénio. Este 
trabalho teve em atenção a influência da radiação, 
incorporando no cálculo um modelo de fluxos [11]. 


Coupland e Priddin [12] aplicaram um modelo 
matemático à previsão do escoamento e combustão 
numa câmara de combustão anular de uma turbina 
de gás industrial. A principal inovação reside na 
geração da malha utilizada. Esta é definida em duas 
dimensões, com o auxilio da transformação de 
Schwarz-Cristoffel e depois gerada em 3 dimensões 
por rotação em torno do eixo de simetria. 
Comparação, com resultados experimentais, do 
campo de velocidades com e sem combustão, da 
distribuição de temperaturas no plano de saída e da 
emissão de NO,, revelou uma precisão aceitável. 
Sampath e Ganesan [13] previram numericamente o 
escoamento e a combustão na mesma câmara de 
combustão e nas mesmas condições experimentais 
usadas por Katsuki, Mizutani e Shibuva [14]. O 
modelo de combustão envolve a resolução de 
equações de conservação da fracção mássica de 
combustível, de CO e da fracção de mistura. Para o 
cálculo da radiação usaram um modelo de fluxos e 
incluíram também um modelo para vaporização e 
combustão de gotas de combustíveis líquidos. Os 
resultados revelaram acordo qualitativo com os dados 
experimentais mas foram observadas diferenças 
quantitativas. Estas foram atribuídas ao facto de a 
malha não ser suficientemente refinada, ao modo 
simplista como foi tratada a injecção de combustível e 
a deficiências do modelo de combustão. 


11) Caldeiras 


Robinson [15] apresenta o primeiro estudo numérico 
tri-dimensional baseado no tratamento integral do 
escoamento, combustão e transferência de calor para 
uma caldeira com queimadores tangenciais de 
grandes dimensões do tipo normalmente usado nas 
centrais produtoras de energia eléctrica. Neste 
trabalho foi utilizado o modelo de turbulência a duas 
dimensões, um modelo de combustão de reacção 
infinita SRQS e o modelo de radiação dos fluxos. O 
algoritmo de cálculo foi só aplicado ao caso de queima 
de combustíveis gasosos. 

Previsões numéricas do escoamento, temperatura e 
fluxos de calor para uma caldeira com queimadores 
tangenciais para combustível gasoso foi apresentado 
por Abbas et al. [16]. O modelo apresentado é 
semelhante ao modelo de [15] sendo no entanto usado 
para a modelação da radiação um sub-modelo muito 
mais poderoso do que o modelo dos fluxos - o modelo 
da transferência discreta - desenvolvido por 
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Lockwood e Shah [17]. Neste trabalho é também 
apresentada a previsão da trajectória de partículas 
não sendo no entanto considerada a sua combustão. 
Tratou-se de um passo importante para o 
desenvolvimento de modelos que considerassem a 
queima de combustíveis sólidos. 

Adbel-Rahman et al. [18] apresentaram um modelo 
matemático para o cálculo do escoamento, 
temperatura, espécies quimicas e fluxos de calor 
numa caldeira a queimar combustível líquido. O 
modelo incorporava o standard modelo de 
turbulência a duas equações, um modelo de 
combustão de reacção infinita e que considerava a 
vaporização instantânea do combustível líquido e o 
modelo da radiação da transferência discreta. O 
modelo foi validado através de comparação com 
dados experimentais obtidas numa caldeira 
semi-industrial. 

Modelos tridimensionais para a previsão de 
escoamento, combustão e transmissão de calor por 
radiação de uma caldeira a queimar carvão 
pulverizado foram apresentados por Boyd e Kent 
[19], Górner e Zinser [20] e mais recentemente por 
Fiveland e Wessel [211. 


Hi) Fornos Industriais 


Os primeiros passos para a formulação de um modelo 
computacional para a previsão dos fenómenos 
ocorrentes num forno tri-dimensional foram dados 
por Patankar e Spalding [1,2]. Neste trabalho 
preliminar, não foi incorporado um modelo de 
turbulência sofisticado, sendo a viscosidade efectiva 
calculada através de uma relação algébrica. 

Modelos físicos mais completos foram incorporados 
no trabalho de Abou Ellail et al. [22] em que foram 
considerados dois fornos semi-industriais, um do Gaz 
de France em Toulouse e outro da IFRF em Ijmuiden, 
Holanda. Neste trabalho foi utilizado um modelo de 
turbulência a duas equações, o modelo de combustão 
SRQS e o modelo de radiação baseado no método dos 
fluxos. 

O trabalho de Pai et al. |23] descreve a aplicação de 


um modelo tri-dimensional no qual foi usado um: 


modelo de turbulência a duas equações, o modelo de 
combustão SRQS e o modelo de radiação dos fluxos a 
um dos fornos semi-industriais da IFRF e a validação 
do modelo através da comparação com resultados 
experimentais obtidos no forno. 

Gosman et al. [24] apresentaram um modelo de 
cálculo tri-dimensional para câmaras de combustão 
de um forno industrial, no qual resolviam o 


escoamento turbulento com combustão e a. 


transferência total de calor, usando para as previsões 
de radiação o método da transferência discreta. As 
previsões estavam de acordo com os pouquíssimos 
resultados disponíveis - temperaturas de saída dos 
gases e temperaturas do tecto. 


Maria da Graça Carvalho 


Carvalho [25] apresentou um modelo para a previsão 
do escoamento tri-dimensional e dos processos de 
combustão e transferência de calor dentro de um 
forno industrial. Neste trabalho foi usado um modelo 
de turbulência a duas equações, um modelo de 
combustão de reacção infinita, um modelo estatístico 
para as flutuações da fracção de mistura, um modelo 
para a formação/oxidação de fuligem e um modelo de 
radiação. O modelo foi aplicado a um forno de fusão 
de vidro com uma geometria complexa e foi 
considerado o caso da queima de vários combustíveis 
gasosos e um combustível liquido. 

Após este trabalho, investigação fundamental 
conducente ao aperfeiçoamento dos sub-modelos com 
especial ênfase no modelo de radiação e de formação 
de fuligem e à inclusão do modelo de formação de NO, 
tem sido realizada no Departamento de Engenharia 
Mecânica do Instituto Superior Técnico. 


1.3 Objectivo da Presente Comunicação 


O objectivo principal da presente comunicação é 
descrever os princípios da modelação física, 
matemática e numérica de câmaras de combustão, 
apresentando cinco exemplos de aplicação a casos 
reais, evidenciando as possibilidades e vantagens da 
aplicação destes modelos no projecto e na procura das 
condições óptimas de funcionamento em câmaras de 
combustão industriais. 


2. DESCRIÇÃO DO MODELO 
FÍSICO-MATEMÁTICO 


2.1 Modelo Físico-Matemático 


O escoamento no interior das câmaras de combustão 
é tri-dimensional, turbulento e, consequentemente, 
não estacionário. Teoricamente seria possível 
resolver as equações exactas, dependentes do tempo. 
Contudo, não é viável a resolução dessas equações 
devido às pequenas escalas de tempo e de 
comprimento associadas à turbulência. Mas do ponto 
de vista da engenharia estamos geralmente 
interessados no conhecimento de valores médios 
temporais das propriedades e não nos detalhes da 
turbulência, pelo que a solução consiste em resolver 
as equações em termos de valores médios. 

As equações de conservação em coordenadas 
curvilíneas ortogonais, podem todas escrever-se na 
forma geral: 
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e 3 dp O íta À 
SPV DO a ENO + S, (1) 


em que & representa a variável dependente. 
Consoante o seu valor assim teremos a equação da 
continuidade, as equações de conservação da 
quantidade de movimento, a equação de 
conservação de energia e a equação de conservação 
de uma espécie química. RE Se 

O tratamento das correlações u'b', v'b' e w'b' exige o 
recurso a um modelo de turbulência. Por sua vez, a 
complexidade do termo fonte, Sg, nas equações de 
conservação de uma espécie química e da energia 
requer o uso de um modelo de combustão e de um 
modelo de radiação que se descrevem seguidamente. 


Modelo de Turbulência 


O modelo de turbulência usado na generalidade dos 
cálculos tri-dimensionais é o k-z [4], baseado na 
hipótese de viscosidade turbulenta isotrópica. 
Envolve a resolução de duas equações de transporte, 
uma para a energia cinética turbulenta, k, e outra 
para a sua taxa de dissipação, £. À análise 
dimensional permite calcular a viscosidade 
turbulenta através de 


2 
C,pk 


= (2) 


E 


em que C, é uma constante do modelo. 
Modelo de Combustão 


O modelo de combustão baseia-se na hipótese de as 
taxas de reacção química serem muito elevadas face 
às taxas de difusão, pelo que prevalece o equilíbrio 
químico. O estado termodinâmico instantâneo do gás 
é descrito em termos de uma propriedade conservada, 
a fracção de mistura, f. Para tal é resolvida a equação 
de transporte para o seu valor médio. 

Em modelação de fornos, a relação entre f e os 
valores instantâneos das propriedades é determinada 
assumindo que a reacção química origina um único 
produto - modelo de reacção química simplificada 
(SRQS) [2]. Em câmaras de combustão de turbinas de 
gás é preferível determinar essa relação com base na 
minimização da energia livre de Gibbs - modelo de 
equilibrio químico completo, Gordon e McBride [71]. 
Em qualquer dos casos, as flutuações da fracção de 
mistura têm de ser consideradas no cálculo dos 
valores médios das propriedades. Desprezar estas 
flutuações resulta na obtenção de perfis de 
temperatura com picos não realistas. Para modelar 
as flutuações é necesséria uma hipótese relativa à 
forma da função densidade de probabilidade (f.d.p.) 
de f (assumiu-se uma distribuição de Gauss truncada 


(26])). A fd.p. de f é especificada em termos do seu 
valor médio e da sua variância - E. 
Consequentemente, é também necessário resolver a 
equação de transporte de g, Spalding [27]. 


Modelo de Radiação 


O modelo de radiação utilizado é o modelo de 
transferência discreta, Lockwood e Shah [17]. 
Baseia-se na solução da equação de transferência de 
calor por radiação para uma dada direcção n 


Ko Tº 
B 


dl 
— = — KI + 
ds n 


(3) 


em que | designa a intensidade de radiação na 
direcção n, T, a temperatura do gás, K o coeficiente 
de absorção do gás, s é a distância percorrida pelo 
raio na direcção considerada e o é a constante de 
Stefan-Boltzmann. 

A equação (3) é resolvida ao longo de raios com 
direcções "representativas" pré-estabelicidas. O 
percurso de cada raio é acompanhado e a intensidade 
de radiação associada ao raio é calculada à entrada e 
saída de cada elemento de volume que o raio 
atravessa. São traçados raios de cada uma das células 
das paredes. À fonte ou poço devido à radiação, em 
cada elemento de volume, é requerida pela equação 
da entalpia que traduz a conservação de energia. 
Para um raio caracterizado por um ângulo sólido 
elementar dt? e atravessando uma área dÃ do 
elemento de volume, perpendicular à direcção 9) do 
raio, a fonte/poço de radiação é calculada através de: 


Srad = (In+1-ln) Q dO dA (4) 


em que 1, e |,,1 São respectivamente as 
intensidades de radiação à entrada e saída do 
elemento de volume. O termo fonte da equação da 
entalpia, para cada elemento de volume, é calculado 
pelo somatório das fontes/poços devidas a todos os 
raios que atravessem esse elemento de volume. 

A emissividade da mistura de gases/fuligem pode ser 
obtida através do modelo de mistura de dois gases 
cinzentos e um gás transparente [28]. 


Modelos de Poluentes 


A combustão de hidrocarbonetos no seio do ar, a 
elevadas temperaturas, origina a formação de 
poluentes, entre os quais NO, e fuligem. À sua 
presença indesejável e o facto da fuligem aumentar 
significativamente a transferência de calor por 


radiação, mesmo quando presente em pequenas. 


proporções, traduz a necessidade de conhecer a 
distribuição espacial das respectivas concentrações. 
Os modelos da fuligem e NO, são descritos de 
seguida. 
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a) Formação/oxidação de fuligem 


A distribuição espacial da média temporal da fracção 
mássica de fuligem - m, - é obtida a partir da 
resolução da sua equação de transporte, na qual são 
introduzidos como fontes, termos referentes à sua 
formação e à sua combustão: S,= Sr-Sa. 

Sr corresponde à taxa de formação de fuligem e é dado 
por Khan e Greeves [29]: 


Sr = Crpru PP exp (- E/Ro T) (5) 


em que pf, é a pressão parcial do combustível, T a 
temperatura dos gases, 4 a razão de equivalência, R, 
a constante universal dos gases perfeitos e CE e n 
constantes do modelo. 

Sa corresponde à taxa de oxidação de fuligem e é dado 
por Magnussen e Hjertager [30]: 


, m 

3. =min[Am (-),A( 

d Ss k ms +m 
5 5 fu 

em que s; e s representam a quantidade de oxigénio 

necessária à oxidação de 1 kg de fuligem e 1 kg de 

combustível, respectivamente, e À é uma constante 


do modelo. 
b) Formação de Nox 


Em sistemas de combustão, com níveis elevados de 
temperatura a formação de óxidos de azoto 
atmosférico de origem térmica é descrita pelo 
mecanismo de Zeldovich: 


N2+ O 


T—* NON 


N+O0 + No+0 


= NO+H 


N + OH (7) 


O0H+0 "“"TTTO+H 


———— 
E 0O+0+M 


M + Os 


A taxa instantânea de formação de NO é dada por 
[31]: 


diNO) 2I0MK K IO JINJ-K K, [NOP) 
dt K, IO, + K [NO] 
em que K representa a taxa de reacção e é obtida de 


Baulch et al. [32]. A resolução da equação de 
transporte da concentração mássica de óxido de azoto 


Ux E 
mem (— (6 
Ra ) 


(8) 
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requer como termo fonte o valor médio da sua taxa de 
formação obtido de: 


“Aa 
Ryo = | Ryo td P (f df (9) 
0 


2.2 Equações do Modelo Matemático 


O modelo matemático descrito nas secções anteriores 
envolve a resolução das seguintes equações 
diferenciais às derivadas parciais: 


- continuidade; 

- conservação da quantidade de movimento: 3 

equações, uma para cada uma das componentes 

da velocidade -U,V e W; 

- energia cinética turbulenta, k; 

- taxa de dissipação da energia cinética 

turbulenta, e; 

- fracção de mistura f; 

- variância da fracção de mistura, Eg; 

- entalpia, h; 
Num sistema de coordenadas curvilíneas ortogonais, 
as equações podem exprimir-se através da forma 
geral [1] em que o coeficiente de difusão To,ef e O 
termo de fonte Sy são discriminadas na Tabela 1. 
As constantes do modelo de turbulência e do modelo 
de combustão, bem como os valores dos números de 
Prandtl turbulentos para as diversas equações, são 
apresentados nas Tabelas ll a IV. 
A transferência de calor por radiação é descrita 
através da equação diferencial ordinária [3]. As 
constantes empíricas necessárias para o cálculo da 
emissividade são apresentadas na Tabela V. 


TABELA I - Termos fonte das equações de 
conservação em coordenadas curvilíneas ortogonais 


p E oret S 

Eq. l 0 

Cont. 

U H,, 
a Lo(V.V) 
a + PE, 3 Per A V)|- 
uv à, v2dl 

Lo PL 
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+ : I (º + E Y + 

| a 1. ay f Hef/0f 0 

1d av vd, vi 

K Y x E Boç'0, gl H, ax ne ày + % e 
3 K, de EA : 

az Bor ax Pur N | dy be, Pp pê 

v al E 1 >( aV ua, ) 

, dx l, dy Per o dx * dy o Ox h HefOh Srad 


TABELA II - Constantes do modelo de turbulência 


d |2 
vas Ar am DO AM] | 
dy » loyla ex CONSTANTE VALOR 
Uº, 15 av  2Uê, ||] 
ret oia (DD) Cu 0.09 
dy Lópl* ef gy , dx C1 1.44 
“al C, 1.92 
E Uvêy 19 (— Nº Ok 0.90 
j | àx 1x yPer ay | ox NA 1.22 
y y y 
U al, à aW' TABELA II - Constantes do modelo de combustão 
“la )| ra (a =| i CONSTANTE VALOR 
que dl, VS V dl, Cgl 
ua, ] dl. “ou va, TABELA IV - Números de Prandtl turbulento para 
- == | 2u (D+) fgeh. 
» ef lay tax ay. 
CONSTANTE 
W u OF 0.90 
id og 0.90 
| 0, 
E ed mms Ha É e e 
"a PE, a ar! V) "E Hof e 
y 
E 1 q l 9V ( aW 
Lloyd ta) my a 
k Hef/0k G-pe 
GQ 0 gs 
E H,/0, 1 - Co prt 
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TABELA V - Constantes do modelo de mistura de 
gases e fuligem. 


n nº Bin B2,nn' x103 kgs ks! 
Óleo pw/pe =1 
l 1 0.717 -0.2964 0.0 925 
1 2 -0.231 0.3861 0.0 4710 
2 l 0.459 -0.1787 2.50 925 
2 2 -0.078 0.1391 2.50 4710 
3 | 0.120 -(0.0499 109.0 925 
3 2 0.0130 -0.0002 109.0 4710 
Gás pw /pc=2 
l 1 0.632 -0.2619 00 1130 
1 2 -0.195 0.3332 0.0 5720 
2 l 0.488 -0.1896 1.88 1130 
2 2 -0.098 0.1844 1.88 5720 
3 l 0.176 -0.0735 68.8 1130 
3 2 0.003 -0.0074 68.8 5720 
Unidades: T IK] 


ken  [m.atmjl 
ksn'  [m.kg/m3)1 


Em — > Qm,n,n' (T) É [1 a exp (- kg, né ks ns) L| 
onde m, é a fracção mássica de fuligem e 
Gm nn' (TD) e Bina + Bo nn CD) 


2.3 Modelo Numérico 


As equações às derivadas parciais constitutivas do 
modelo matemático são discretizadas pelo método do 
volume de controlo conduzindo a equações às 
diferenças finitas. A discretização do termo 
convectivo é feita através do esquema híbrido |33]. 
Para resolver essas equações subdividiu-se o domínio 
de cálculo num número finito de volumes de controlo 
- as células. As velocidades e as pressões são 
calculadas através de uma variante do algoritmo 
SIMPLE [34] para o caso dos cálculos com combustão, 
sendo usado para os cálculos sem combustão 
(exemplo, tanque de vidro em fusão) o algoritmo 
PISO [35]. 
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A solução dos sistemas de equações algébricas 
individuais é obtida combinando o algoritmo TDMA 
("Tri-diagonal matrix algorithm") com o método de 
Gauss-Seidel. 


3. RESULTADOS E DISCUSSÃO 


A aplicação dos modelos atrás descritos a diversas 
geometrias em condições de funcionamento reais, é 
descrita de seguida. 


3.1 Câmara de Combustão de uma Turbina de 
Gás. 


Nas duas últimas décadas, a consciência da escassez 
das fontes de energia, acompanhada do aumento do 
preço dos combustíveis, por um lado, e a legislação 
limitando a emissão de poluentes libertados na 
combustão, por outro, conduziram à necessidade 
permanente de optimizar o funcionamento dos 
equipamentos onde se dá a combustão, entre as quais 
as câmaras de combustão das turbinas de gás. 

Face à multiplicidade dos requisitos a satisfazer, 
interdependentes e alguns deles incompatíveis 
(eficiência de combustão elevada, ignição fiável, 
limites de flamabilidade amplos, insensibilidade a 
manifestações de instabilidade na combustão, baixa 
queda de pressão, distribuição de temperaturas à 
saída projectada de modo a maximizar a vida das pás 
da turbina, custo mínimo, durabilidade máxima, 
facilidade de manutenção e baixas emissões de fumo 
e de poluentes) a optimização do funcionamento só é 
possível mediante um conhecimento detalhado dos 
fenómenos que ocorrem no interior da câmara de 
combustão. Daí o esforço que tem sido desenvolvido 
nos últimos anos no sentido de melhorar esse 
conhecimento, tanto a nível experimental como a 
nível da modelação matemática. 

A despeito do numeroso trabalho desenvolvido na 
modelação matemática de câmaras de combustão de 
turbinas de gás, referido na secção 1.2., tem sido 
prestada pouca atenção à transmissão de calor por 
radiação. A maioria dos autores que usaram modelos 
matemáticos para o estudo de turbinas de gás 
ignoraram a radiação, com o argumento de ser 
desprezável face ao calor libertado na combustão. No 
entanto, se quisermos conhecer a distribuição de 
temperaturas nas paredes da câmara de combustão, o 
cálculo dos fluxos de radiação é indespensável. O 
interesse prático desse conhécimento advém da 
circunstância de as paredes não poderem suportar, 
para os materiais de uso corrente, temperaturas 
superiores a 1100 K [36]. Dado que na zona primária 
a temperatura do gás pode atingir valores da ordem 
de 2300 K é necessário que as paredes sejam 
arrefecidas, de modo a garantir que não se ultrapassa 
aquele limite. Caso contrário, as tensões de origem 
térmica poderão levar à propagação de fissuras e 
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consequente redução do tempo de vida do "liner" 
(invólucro da câmara de combustão). 

Por outro lado, um cálculo rigoroso da emissão de 
poluentes necessita também de um conhecimento 
preciso da transferência de calor por radiação, uma 
vez que os modelos para e cálculo da emissão do 
monóxido de carbono, óxidos de azoto ou 
hidrocarbonetos não queimados requerem um 
conhecimento preciso da distribuição de 
temperaturas em toda a câmara. 

Além disso, nos últimos 30 anos tem-se assistido a 
um aumento da razão de compressão nas turbinas de 
gás e parece provávi| que essa tendência se 
mantenha. Isto é acompanhado por um aumento de 
calor transferido por radiação, justificando a 
importância de que o fenómeno se reveste. Tornou-se 
portanto premente desenvolver modelos de radiação 
ou adaptar os modelos já existentes às condições das 
câmaras de combustão de turbinas de gás. 

O conjunto dos modelos matemático e fisicos de 
turbulência, combustão (equilibrio químico) e 
radiação atrás descritos, foram englobados num 
algoritmo de cálculo numérico que foi aplicado a uma 
câmara de combustão anular com 18 queimadores 
(Figura 1), em que o combustível é propano. As 
paredes anulares interior e exterior têm orifícios 
para a entrada de ar secundário e de diluição. Na 
secção de entrada, junto às paredes anulares, é 
admitido ar de arrefecimento que se destina a 
impedir que a temperatura dessas paredes atinja um 
nível tal que as tensões de origem térmica 
provoquem a danificação permatura da câmara. 


Figura 1 - Geometria da câmara de combustão da turbina de gas. 


No trabalho de Carvalho, Durão e Lockwood [37] foi 
calculada a distribuição dos fluxos por radiação 
através do uso do modelo de transferência discreta. A 
equação da entalpia foi considerada e nela foi 


incluido um termo de fonte devido à radiação. À 
distribuição de temperatura foi calculada a partir da 
distribuição de entalpias. No entanto, neste trabalho, 
a distribuição de temperaturas nas paredes foi 
assumida (tomada como constante e igual à 
temperatura do ar de admissão), a qual foi usada 
como condição de fronteira para o cálculo da 
transferência de calor por radiação. 

O trabalho de Carvalho e Coelho [38] dá um 
contributo importante para a modelação dos 
fenómenos ocorrentes dentro das câmaras de 
combustão de turbinas a pás, ao apresentar pela 
primeira vez na literatura o cálculo simultâneo da 
distribuição dos fluxos e das temperaturas nas 
paredes. Para conseguir determinar a distribuição de 
temperaturas e de fluxos nas paredes da câmara, foi 
estabelecido um balanço energético aplicado a cada 
elemento de volume adjacente à fronteira do domínio 
de cálculo. 

No trabalho de Carvalho e Coelho [39] foi investigado 
o efeito da pressão na transferência de calor por 
radiação para três condições de funcionamento 
(Tabela V]). 

Na figura 2 são apresentadas, como exemplo, as 
distribuições de fluxos e de temperaturas na perede 
anular superior para a pressão de 25 bar. 


TEMP. NAS PAREDES [K) 
PAREDE ANULAR EXTERIOR 


FLUXO RADIAÇÃO (KW/m?) 
PAREDE ANULAR EXTERIOR 


Figura 2 -Distribuição de fluxos e de temperaturas da parede 
anular superior (p = 25 bar) na câmara de combustão da turbina 
de gás. 


A localização dos orifícios de entrada de ar 
secundário e do ar de diluição é facilmente 
reconhecível pela presença de isotérmicas fechadas 
de baixas temperaturas. Na zona primária e no início 
da zona intermédia, a temperatura é bastante 
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Pressão (Bar) 


Temp.entrada 
ar (K) 


Temp.entrada 
combustível (K) 


Caudal de ar 
primário (Kg/s) 


Caudalar 
arrefecimento 


(Kg/s) 


Caudal ar secun- | 


dário (Kg/s) 


Caudal ar dilui- 


ção (Kpg/s) 


Caudal combus- 
tível (Kg/s) 


Razão 
ar/combustivel 


TABELA VI - Condições de funcionamento da 
câmara de combustão da turbina de gás. 


uniforme, inferior a 1000K, evidenciando o efeito 
protector dos jactos de ar de arrefecimento. À medida 
que se caminha para jusante, a difusão provoca a 
atenuação progressiva desse efeito protector, pelo 
que a temperatura na perede tende a aumentar. 

O conhecimento da distribuição de temperaturas nas 
paredes reveste-se de grande interesse prático 
porque, como já foi referido, a temperatura das 
paredes não deverá ultrapassar um determinado 
nível, que para os materiais de uso corrente é 1100K. 
Verifica-se que no caso estudado, esta restrição é 
razoavelmente satisfeita, com excepção das regiões 
da zona de diluição onde a temperatura excede os 
1150K. 

Os fluxos de radiação aumentam progressivamente 
na zona primária para jusante, desde valores da 
ordem dos 100 KW/m2 até fluxos da ordem dos 200 
KW/m2. Na zona intermédia os fluxos atingem 
valores elevados, chegando em alguns locais a 
ultrapassar os 300 KW/m>2, o que é consequência das 
temperaturas elevadas dos gases nesta região. Na 
zona de diluição a tendência é no sentido de uma 
diminuição dos fluxos devida ao elevado valor das 
temperaturas da parede. 
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Nota-se também a ocorrência dos valores mais 
elevados dos fluxos nos locais dos orificios de ar 
secundário e de diluição. Estes fluxos são mais 
elevados do que para as paredes mas, dado que os 
orifícios cobrem uma área inferior a 5% da área total 
da parede, este facto não é significativo. 

Na Tabela VII são apresentados os fluxos de 
radiação, convecção e total nas paredes da câmara de 
combustão para os 3 casos estudados, bem como os 
quocientes entre esses fluxos e a potência calorífica 
libertada na combustão. Pode verificar-se que os 
fluxos são pequenos em face da energia libertada na 
combustão. Os fluxos de convecção são ligeiramente 
superiores aos de radiação mas, graças aos jactos de 
ar de arrefecimento, mantém-se em níveis baixos. 
Este ar de arrefecimento reduz a eficiência da 
combustão e contribui para a presença de poluentes 
nos gases de escape mas a sua presença é necessária 
para manter a temperatura das peredes num nível 
aceitável. 


Qr fluxo de calor nas paredes devido à 
radiação. 

Qe fluxo de calor nas paredes devido à 
convecção. 

Qt+Qr=Qe fluxo total nas paredes. 

P potência calorífica libertada na 


combustão. 


Neste estudo, a razão ar/ combustivel foi mantida 
constante, sendo aumentada a pressão, a 
temperatura e o caudal do ar de admissão. 

O aumento da pressão traduz-se num aumento 
substancial dos fluxos, que passaram para cerca do 
dobro quando a pressão passou de 5 a 15 bar, Mas o 
aumento de 15 para 25 bar levou a um aumento mais 
ligeiro dos fluxos. Esse aumento á cerca de 10% para 
a radiação e ainda menor para a convecção. 

Os resultados destes trabalhos [37,38,39] 
demonstraram que a distribuição de temperaturas 
nas paredes pode ser calculada, para qualquer 
câmara, através de um balanço energético às 
paredes. O método delineado constitui o único 
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processo correcto para o cálculo dessa distribuição de 
temperaturas. 


3.2 Forno de Vidro de Garrafas da Sociedade 
Saint-Gobain 


Os modelos anteriormente descritos, à excepção do 
modelo e formação de óxido de azoto, foram aplicados 
a um forno industrial de fusão de vidro para garrafas 
da Sociedade Saint-Gobain, [40,41,42,43). 

Este forno (vide Figura 3) é aquecido com uma chama 
em ferradura com dois pórticos colocados numa das 
paredes. Cada um dos pórticos contém 5 injectores. O 
forno queima alternadamente a partir de cada um 
dos pórticos, funcionando cada pórtico de entrada ou 
de exaustão. Os produtos de combustão passam 
através de um regenerador que serve para 
pré-aquecer o ar da combustão antes dele entrar no 
forno. O ar é injectado com uma inclinação de 16º em 
direcção à soleira do forno.e o combustível com uma 
inclinação de 8º em direcção ao tecto do forno, sendo 
estes valores ajustáveis. As matérias primas (areia, 
aragonite, dolomite) e o vidro reciclado entram no 
forno, na parede lateral junto à parede que contém os 
pórticos e flutuam à superfície do vidro fundido. Esta 
carga sofre reacções químicas, funde e vai ser 
misturada por convecção natural e difusão. O vidro é 
removido através de uma garganta existente na 
parede oposta à parede que contém os pórticos. 


O principal objectivo deste estudo foi determinar a 
geometria óptima do forno e as condições óptimas de 
funcionamento. O referido estudo permitiu a análise 
das condições de funcionamento do forno a queimar 
gás natural e fuelóleo pesado. 


A simulação foi dividida em três partes, sendo a 
primeira referente à simulação da câmara de 
combustão (para gás e para óleo). Para a simulação 


= E s  —— — = "— — em 


da região de fusão da matéria prima, foi usado um 
programa bi-dimensional parabólico com um modelo 
para a reacção e outro para a fusão. Por fim, para a 
simulação do tanque com vidro fundido foi usado um 
programa elíptico laminar tri-dimensional. 

A figura 4 mostra as distribuições do fluxo de calor 
para o vidro nos casos em que o combustível é um 
fuelóleo e gás natural. Verifica-se que queimando 
fuelóleo o fluxo de calor para o vidro é superior. Às 
diferenças verificadas na performance do forno 
queimando gás natural e fuelóleo (Tabela VIII), são 
devidas à maior emissividade da chama no caso do 
fuelóleo tendo como consequência o aumento de 
transferência de calor para o vidro. Os fornos 
queimando gás natural têm maior duração devido à 
temperatura da abóbada ser menos elevada e mais 
uniforme. Os fornos queimando gás natural são 
limpos e conduzem a uma menor emissão de 
poluentes. 


Os resultados obtidos concordam (qualitativamente) 
com os poucos resultados experimentais existentes. 


3.3 Forno de Vidro Cerâmico da Metal 
Portuguesa 


A presente aplicação consiste na simulação numérica 
dos fenómenos ocorrentes na câmara de combustão de 
um forno de vidro cerâmico, utilizando os modelos 
atrás descritos, excepto o modelo de emissão de NO,, 
[44, 45]. 

Com o fim de validar os métodos de cálculo e 
aprofundar o conhecimento dos processos de 
combustão foi desenvolvido um programa de 
investigação experimental no forno em condições 
normais de funcionamento | 46, 47]. 

Conforme a Figura 5, trata-se de um forno de 
funcionamento contínuo com escoamento em 
ferradura, em que o comburente é o oxigénio e o 


1- CÂMARA DE COMBUSTÃO 

2 - TANQUE DE VIDRO 

3- PÓRTICO DE ENTRADA 

4 - PÓRTICO DE SAÍDA 

5 - QUEIMADOR COM 5 INJECTORES 

6 - DEGRAU 

7- GARGANTA 

8- ALIMENTAÇÃO DE MATÉRIA-PRIMA 


Figura 3 - Geometria do forno Saint-Gobain. 
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A=221520 
B=189850 
C =158180 
D=126510 
E =94840 


F = 63170 
6 = 31490 
H=-170 

| =-31840 


gds natural 


A=445010 
B = 309060 
C = 260520 
D = 186180 
E = 123770 


F=73300 
6=34770 
H= 10900 
| =-=32420 


fuel óleo 


Figura 4 - Distribuição do fluxo de calor para o vidro em [W/m?] 
no forno de Saint-Gobain 


combustível um fuelóleo pesado. O oxigénio é 
admitido no forno à temperatura ambiente enquanto 
que o combustível é pré-aquecido. Existe apenas um 
queimador descentrado que se encontra inclinado 
para o chão a 4º. A porta de saída dos gases de 
combustão, que comunica directamente com a 
chaminé de exaustão dos mesmos, serve de saida 
para a matéria em fusão e encontra-se situada numa 
parede lateral. O forno é alimentado de matéria 
prima por uma porta situada na parede oposta à do 
queimador, originando aí um empilhamento desse 
material por acumulação do mesmo. O escoamento 
da matéria fundida para o exterior, faz-se por efeito 
da gravidade, em contra-corrente com a chama. 

As condições de funcionamento do forno 
encontram-se na Tabela IX. 
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Gás i 
Estaca Fuelóleo 
Fluxo de calor por 

radiação para a 
superfície de vidro 1.0661 x 
[W/m2] | 107 


1678 


Fluxo de calor por 
convecção para a 
superfície de vidro 
[W/m2] 


Temperatura dos 
gases à saída 


Temperatura máxima 


da abóboda [K] 


Velocidade de 
saída/m/s | 


forno da Saint-Gobain queimando gás natural e 
fuelóleo. 


EM 
Po 
* 1 - QUEIMADOR 


e 

2 - PÓRTICO DE EXAUSTÃO 
3 - ALIMENTAÇÃO DE MATERIAS-PRIMAS 
4- SAÍDA DE CARGA 


gt 
e 
| . 


YU) 
EV 


Figura 5 - Geometria do forno da Metal Portuguesa. 


A figura 6 mostra, a título de exemplo, a comparação 
entre os perfis de temperatura e as concentrações de 

oxigénio medidos e calculados. Como se pode 
verificar, a concordância entre os valores é excelente. 
À temperatura média do gás e a concentração de 
fuligem estão representadas na Figura 7 para um 
plano vertical que atravessa o queimador. À 
temperatura média de gás é praticamente 
homogénea numa grande extensão do forno (valores 
variando entre 1600K e 1700K). Na zona do 
queimador a temperatura atinge 3000K. Tomando 
em consideração que a temperatura adiabática da 
chama para as condições oxifuel excede 8000K, estes 
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Figura 6 - Comparação entre valores de temperaturas dos gases e 


concentração de oxigénio medidos e calculados no forno da 
Metal Portuguesa. 


1.3 1.8 


resultados mostram a dominância da transferência 
de calor por radiação. 
De igual modo as concentrações mássicas de fuligem 


são praticamente homogéneas e com um valor perto 


de zero em grande extensão do forno. 


6 Caudal de combustível = 0.03 Kg/s 
6 Caudal de oxigénio = 0.105 Kg/s 


e Excesso de oxigénio = 2% 
e Temperatura de entrada do combustível = 393K 
e Temperatura de entrada do oxigénio = 300K 


"TABELA IX - Condições de funcionamento do 
forno da Metal Portuguesa 


Figura 7- Campos escalares num plano vertical que atravessa o 
queimador no forno da Metal Portuguesa. 

a) Temperatura média dos gases 

b) Concentração mássica de fuligem 


3.4 Forno de Vidro de Placa Plana Tipo 
"Pittsburgh" 


Este estudo apresenta a simulação completa de um 
forno de vidro industrial. O cálculo da câmara de 
combustão é feito de um modo semelhante ao cálculo 
efectuado no caso dos fornos da Saint-Gobain e da 
Metal Portuguesa, no entanto um novo modelo 
bastante mais elaborado é usado para o cálculo dos 
fenómenos ocorrentes no tanque de vidro [51, 52). 
Trata-se de um forno de vidro placa plana do tipo 
Pittsburgh com queimadores laterais (Figura 8). O 
forno tem 4 pórticos em cada uma das paredes 
laterais e cada pórtico tem dois queimadores. O forno 
queima alternadamente de cada lado, com intervalos 
de uma hora. Os produtos de combustão passam 
através de regeneradores sendo o ar de combustão 
pré-aquecido. O forno queima um fuelóleo pesado e o 
comburente é ar. 


As condições de funcionamento encontram-se na 
Tabela X. | 
Os principais objectivos deste trabalho foram a 
investigação das correntes de convecção e a 
localização do "ponto-quente" no tanque de vidro em 
fusão e, a procura de uma melhor eficiência de 
combustão, uma vez que se verificou que 0 processo 
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MATERIAL EM FUSÃO 


QUEIMADORES 
/ PÓRTICO DE ENTRADA / SAÍDA 


CÂMARA DE COMBUSTÃO —T | 
GARGANTA 


TANQUE DE VIDRO — 


Figura 8 - Geometria do forno de vidro de placa plana tipo 
“Pittsburgh”. 


PÓRTICO N Nº 
Caudal total de | 
combustível ei 496 N to | 
Caudal total dear | | | 
[Kg/h] 5480 | 8100 | 8680 | 


emaum [u[n [o [8 


Velocidade do 
0. 505 | 0. 765 O. os 
Velocidade doar | | 

[M/s] | 


TABELA X - Condições de Ee do forno 
de vidro de placa plana tipo "Pittsburgh" 


combustivel [M/s] 


de mistura ar/combustível era pobre, existindo ainda 
combustível não queimado à saída do forno. 

A figura 9 contém a distribuição de temperaturas no 
tanque de vidro mostrando que os valores elevados 
da temperatura e dos seus gradientes ocorrem na 
região do "ponto-quente” (com um valor máximo de 
1630 K) que se localiza próximo da região do terceiro 
pórtico de entrada. 

Como atrás se referiu, foi efectuado um estudo 
paramétrico com o intuito de melhorar as condições 
de funcionamento da câmara de combustão [49, 50). 
Assim, estudou-se o efeito do pré-aquecimento e do 
excesso de ar e o efeito da posição do queimador. 
Apenas este último revelou melhorias significativas 
na eficiência do forno. 


Maria da Graça Carvalho 


A figura 10 mostra a distribuição da fracção de 
mistura num plano vertical, que contém o 
queimador, normal aos pórticos de entrada e 
saida,para o caso real (caso 1) e para o caso em que a 
distância do queimador ao chão do forno foi duplicada 
(caso 6). Verificou-se que o processo de mistura dos 
reagentes melhorou na medida em que deixou de 
haver combustivel à saida do forno e o oxigénio não 
consumido, anteriormente existente na região 
superior da câmara de combustão, passou a interferir 
no processo de mistura. 

Os resultados obtidos concordam com os poucos 
valores experimentais para este forno. 

Uma conclusão importante deste trabalho foi o facto 
de os fenómenos ocorrentes na câmara de combustão 
e no tanque de vidro estarem intimamente 
interligados e, a não ser que hajam medições precisas 
da temperatura da superficie de vidro, os cálculos da 
câmara de combustão e do tanque de vidro deverão 
ser feitos em conjunto e interligados. 


Para um forno de vidro de placa plana semelhante ao 
apresentado na Figura 8, foi elaborado um trabalho 
[51] com especial ênfase no problema da emissão de 
poluentes. Os modelos físico e matemático utilizados 
foram os mesmos do que no caso descrito 
anteriormente, tendo no entanto sido incorporado o 
modelo de formação de NO. Por razões de sigilo, as 
condições de funcionamento do forno não são 
apresentadas. Pela mesma razão, os resultados são 
apresentados numa forma adimensional. 

A concentração mássica de NO é mostrada na Figura 
11. Os resultados mostram que o NO térmico é 
produzido principalmente no bordo exterior da 
chama, na região próxima da mistura 
estequiométrica, onde as temperaturas são elevadas 
e a concentração de oxigénio é considerável. A 
concentração de NO é máxima perto da frente de 
chama, diminuindo progressivamente para jusante 
devido à convecção e à difusão para longe da região 
onde é formado. No pórtico de saída a concentração 
mássica de NO é uniforme. Este valor está em bom 
acordo com o valor medido na fábrica. 


3.5 Caldeira 


O conjunto dos modelos físico-matemáticos descritos 
na Secção 2 foram englobados num algoritmo de 
cálculo numérico que foi aplicado a uma caldeira de 
uma central de produção de energia de modo a 
calcular o escoamento e a transferência de calor para 
as paredes da caldeira. 

A caldeira examinada no presente estudo foi 
convertida de queima de carvão para queima de óleo 
e tem um débito de vapor máximo de 75 Kg/s a 66 bar 
e a 490ºC. As principais dimensões da câmara de 
combustão estão indicadas na Figura 12. A caldeira 
tem 16 queimadores, 8 em cada parede lateral. 
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No trabalho de Carvalho e Coelho [52] foi calculado o 
campo de velocidades, temperatura, concentrações 
das diversas espécies químicas e fluxos de calor para 
as paredes da caldeira. 
Os resultados foram comparados com os resultados 
experimentais de [53]. 


A título de exemplo, a Figura 13 apresenta a 
distribuição de fluxo de calor total calculado e 
medido para uma das paredes da caldeira. Os valores 
do fluxo de calor concordam com os valores 
experimentais disponíveis. 


4. CONCLUSÃO 


O presente trabalho descreve um algoritmo de 
cálculo do escoamento, transferência de calor, 
combustão e emissão de poluentes no interior de 
câmaras de combustão reais. O algoritmo de cálculo 
foi aplicado a cinco câmaras de combustão: câmara de 
combustão de uma turbina de gás; um forno de vidro 
para garrafas; um forno de vidro cerâmico; um forno 
de vidro de placa plana e uma caldeira do tipo 
utilizado em centrais produtoras de energia. 

As previsões efectuadas com o algoritmo atrás 
referido, para o forno de vidro cerâmico da Metal 
Portuguesa, foram extensivamente validadas com 
resultados experimentais medidos no forno em 
funcionamento normal, tendo-se verificado uma 
excelente concordância entre eles. 

Este algoritmo constitui um poderoso meio não só 
para o projecto de engenharia de câmaras de 
combustão como também para a procura das 
condições óptimas de funcionamento de sistemas de 
combustão já existentes. 


TEMPERATURA [(K) 


A = 0.80 
B = 0.50 
C=0.2 
D = 0.10 
E = 0.05 
F = 0.025 


Figura 10 - Distribuição da fracção de mistura na câmara de 
combustão (caso 1 e caso 6) do forno de placa plana tipo 
“Pittsburgh”. 
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SUMMARY 


Spectral Analysis comprises three major steps: 


methodological support, algorithmic implementation 


luation procedure 


The first step deals with the theoretical foundations 
of the spectral estimation models. These, as it is 
shown, must provide for objective and coherent use of 
the available information. 

The algorithmic implementation step supplies the 
user with the algorithms to implement the methods 
found in step one. On the other hand, mathematical 
justification and some properties of the algorithms 
must be stated in this step. 

The evaluation procedure step allows users to chose, 
among several spectral estimators, one to use with 
their particular data problems, according to the 
constraints of the data space and the information to 
be obtained from spectral estimates. 

This paper discusses this three step procedure for 
spectral estimation with emphasis on the general 
framework for the evaluation of spectral analysis 
techniques from the point of view of the user. 


1. INTRODUCTION 


Two decades of active research led to a multitude of 
methods and algorithms to perform spectral 
estimation. However, there is no general guide to 
help a user, having a particular data problem, in the 
difficult task of performihg an adequate spectral 
estimation. We can classify user's problems in two 
types: a) selection of valid estimators and b) choice of 
the adequate estimator for his data problem. His 
difficulties in selecting valid estimators come from 
the following facts: 


1- The methods are not univocally labelled. 

2- Some methods are not valid, since they lead 
to non-coherent estimates. 

3- There is no clear distinction between method 
(model) and algorithm. 

4- Some methods have several implementation 
algorithms. 

5 - Some of the algorithms used to implement a 
given method may be equivalent although 


having different formulation and 
mathematical support. 
6 - Some estimation procedures are "ad hoc" and 


do not have a correct mathematical support. 


In this paper we will show how to overcome those 
problems. This is the object of the first and second 
procedure steps for Spectral Analysis described in 
sections 2 and 3, respectively. The first step is called 
mathematical support and supplies a user with 
spectral analysis models. The second step is 


. concerned with the algorithmic implementation and 


allows him to find correct algorithms to implement 
the methods obtained in the first step. As a 
consequence of this formulation we define any 
spectral estimator in terms of a pair 
method-algorithm. 

Supose now that our user followed these steps and 
reached to a situation of having several spectral 
estimators. Thus, he is faced with another problem: 
How to choose one for his particular data problem? 
The answer to this question is the objective of the 
third step of Spectral Analysis: the evaluation. This 
is a very important problem which has not received 
enough attention in the literature. 

Several attempts have been made to pursue it but 
from very restrictive points of view and with no 
useful final results. For example, in connection with 
parametric models several tests have been proposed 
to verify the adequacy of a given model to the signal 
being used. The most important of them is the Q-test 
(Quenouille, 1947) for AR signals. Within the scope 
of a specific spectral analysis method several tests 
were proposed which, directly or indirectly, make 
comparisons between the correct and estimated 
spectra (see [4], and [6]). 

In this paper and following [13], we adress the 
problem of evaluating spectral estimators according 
to some pre-specified criteria and working space of 
signals. In other words, we suppose given a class of 
stochastic processes and a criterion (or criteria) 
under which we want to measure each spectral 
estimator, as estimator for the class. 
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Of course, we do not intend to give an absolute 
answer for this question, but to formulate a general 
framework which allows the possibility of finding 
coherent solutions to the problem. The requirements 
for what we intend to be a good evalution scheme are 
presented in sec.4.1 In sec.4.2 we describe one of the 
frameworks for the evaluation presented in [11]. 
Simulation results illustrating the proposed 
formulation are shown in sec.5. Sec.6 is devoted to 
some conclusions. 


2. METHODOLOGICAL SUPPORT 


This first step for Spectral Analysis supplies answers 
to the following statement: 
Compute the spectrum, S(f), of a stationary 
stochastic process Xn from some "a eita 
: | | Rn Ê 


about the process. 
In order to approach a solution to the problem we 


must realize that |2]: 


a) The spectrum has to be consistent with the 
available information and to match it perfectly. 

b) No “a priori" assumption about nonavailable 
information should be made. That is: the 
spectrum must be non commital in relation to 
unavallable values. 


It seems clear that there are an infinite number of 
solutions to the problem. The selection of one is a 
subjective task. Burg |2] proposed a variational 
principle: 

Find an extremal of a functional of the spectrum 

under certain constraints, set by the "a priori" 

information. 
The choice of the functional is done according to 
other kind of considerations: Mathematical, 
Physical, Philosophical, etc. Sometimes 
combinations of functionals are used in order to 
achieve specific goals (see [17, 18]). Some of the most 
important methods derived from this framework are: 
Maximum Entropy Method [2], Maximum Flatness 
Method [10, 12], Minimum Cross-Entropy Method 
[15] and Minimum Itakura-Saito Distance Method 
[41]. 
It is interesting to remark here that two famous 
spectral estimation methods are not correct 
“vis-a-vis" our approach: conventional and Capon's 
ML methods. Conventional method assumes a zero 
extrapolation of the ACF and Capon's ML method is 
not consistent with the “a priori" information, since 
the inverse Fourier transform of the ML spectrum is 
different from the Autocorrelation Function used as 
starting point. 


3. ALGORITHMIC IMPLEMENTATION 


In the previous section we showed how different 
methods can be constructed from the "a priori” 
information. However, nothing was said about the 
acquisition of that information and its use in the 
computation of spectral estimates. This is the 
problem to be faced in this section reformulated in 
the more realistic statement: 

Estimate the spectrum S(f) of a stationar 

realization. E 

This means that we must estimate the parameters 
needed to define the spectrum estimate from a finite 
realization of the stochastic process. Globally we can 
group the algorithms in two classes according to the 
use of the signal: 


C1 - Estimate the ACF and, from it, the spectral 
parameters. 
C2 - Estimate the spectral parameters by looking 
for an extremum of a functional involving some 
relations between the signal and those 
parameters. 


In the first case, the parameters may be estimated by 
two different procedures: 


a) - Using theoretical relations between the ACF 
and those parameters. 

b)- Making an interative adjustment of the 
parameters, by minimizing a function of the 
error between the original and actual ACF. 


Certain algorithms modify procedure a) to include 
ACF estimation errors and minimize an error 
functional. In this case we must refer the use of the 
SVD algorithms (Cadzow [20]). 

The algorithms of class C2 do not have a so wide 
application as those in class Cl, because they have 
been developed mainly for rational methods. In this 
case, maximum likelihood techniques and 
least-squares algorithms [1], [8-10] are currently 
used. Some of the best known algorithms can be 
easily classified according to those classes, namely: 
The algorithms presented in [10] and [11] to 
implement maximum flatness method belong to class 
Cl: 

The Yule-Walker algorithm for maximum entropy 
AR method belongs to class Cl. 

The well known Burg method, used in AR spectral 
estimation, belong to class C2. 

The algorithms we just referred are used in spectral 
estimation of regular stationary stochastic processes. 
In the case of singular precesses some special 
algorithms are available: Prony method [7], 
Pisarenko harmonic decomposition [7] and 
Kumaresan-Tufts SVD algorithm [17]. 
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4.1 REQUIREMENTS FOR THE EVALUATION 


Consider the space of the stationary stochastic 
processes and suppose it devided into classes 
according to some caracteristics of the processes |e.g. 
the class of the AR(N) processes]. Take out a process 
from one class; which spectral estimator should we 
use to estimate the spectrum from finite realizations 
of the process? 

The methodology for answering the previous 
question considers each spectral estimator as a 
“black box” (having finite realizations as input and 
spectral estimates as output) and lies in four items: 


a) Selection of certain features computable both 
from the spectrum and from the signal. 

b) Algorithms for the computation of such 
features. 

c) A measure of "distance" between features. 

d) A decison rule. 


“The selection of features depends on the chosen 
criterion and may be the most difficult task in the 
evaluation procedure; the answers to the other three 
items arise almost naturally from the caracteristics 
of the features - see the appendix for an example. 

The search for features according to a given criterion 
remains to be the main challenge of the evaluation 
framework. The higher order zero-crossings, 
reported in the appendix, áre suitable for the 
evaluation under bias. It would be important to find 
features according to variance or resolution. 

In the following section, we suppose the four previous 
requirements to be satisfied. 


4.2 THE FRAMEWORK 


In fig. 1, we propose a framework for the evaluation 
of spectral estimators. The procedure for two spectral 
estimators is ilustrated; the generalization is simple. 


The general procedure will be: 
1-Fork=1,...,K 


a) Generate a finite realization of the process. 

b) Extract its features. 

c) Compute the corresponding spectral 
estimates Sm() (m=1,...,M). 

d) From each spectral estimate, Sm(f), compute 
the features. 


2 - Compute the mean values of the features in 1 - b). 
3 - Compute the mean values of the features in 1 - d). 


or, ifthe true spectrum is known, 
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3 - Extract those mean values from the spectrum. 
4- For each m=1,...,M measure the “distance” 
between the values computed in 2 and 3 or 3'. 


We choose the spectral estimator with least 
“distance”. 


FE - feature extractor 
A - average 
D - distance 


Fig.1- SG-signal generator 
SE - spectral estimator 
C - comparator 
TS- true spectrum 


6. SIMULATION RESULTS 


The illustration of the three step procedure we have 
just described will be done in this section using a 
simple example. In the first step we maximize the 
entropy under correlation constraints to obtain 
maximum entropy AR method. It is known to exist 
several algoritms to implement this method [13]; we 
selected four of them: modified covariance (MC) 
algorithm [13], LSQ Marple's algorithm [9],Burg 
technique (BT) [2,3] and a modified Burg technique 
(MBT) [13]. 

Table 1 shows the results obtained from the 
comparison of the four estimators, using bias as 
evaluation criterion, Here, we use the higher order 
zero-crossings [eqgs. (5) and (6) in appendix |] as 
features and the statistic dy, [eq. (8) in appendix] as 
"distance". The use of these features was motivated 
by the knowledge of the one to one relationship 
between the zero-crossings and the ACF, We used 
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100 50-point long realizations of an AR(4) gaussian 
stochastic process. We present the results obtained 
with (first line) and without (second line) the 
knowledge of the correct spectrum. It is important to 
remark the coincidence between these simulation 
results and the theoretical ones presented in [13]. 


7.-CONCLUSIONS 


We presented a three step formulation of Spectral 
Analysis. As shown, this formulation supplies a user 
with a well defined framework which allows him to 
solve his particular data problem. We made a 
distinction between method and algorithm and 
defined spectral estimator. 

Each spectral estimator is the result of the first and 
second steps in the present formulation: 
Methological Support and Algorithmic 
Implementation. These steps put at the disposal of a 
user several estimators. 

A special attention was directed towards the 
quantitative evaluation of spectral estimators from a 
user's point of view. We presented a general 
framework for such evaluation. The proposed scheme 
has two important caracteristics: it is easily 
implementable and does not need knowledge of the 
true spectrum. 

We reported some simulations results concerting the 
application of the procedure. 
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APPENDIX 
THE HIGHER ORDER ZERO-CROSSINGS 


Let z, be a zero-mean stationary stochastic process 
and let D denot the backward difference operator: 


Dz =2 2. (1) 


Define the clipping operator U by: 


U z= | é (2) 
n -1 if 


and 


x = 1) a (4) 


Suppose that z, and x, observed during a L-length 
interval. The numbers Ny defined by: 


L 
> [xi | (5) 


n n-l 
n=1 


IRA 
“e 5 
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are called higher order zero-crossings of z,. TABLE 1 


The expected value of Ny) is not easy to compute in 
the general case. However, if z, is à gaussian process 
(13): 


- L=1 P) = El (6) 
e[N [== iam 


where,j=0 and 


System: AR (4)- Poles: .8 + 3.:.66 + 142. 
R, (n) being the autocorrelation function of zh. As it GAUSSIAN SIGNALS 

can be seen, for a fixed j, E[N,)| depends on the j+2 

initial values of the autocorrelation function. It is not 

hard to conclude that there is a one to one 

relationship between the set of autocorrelation 

values and that of the higer order zero-crossings. 


Asin [11], defines: 


J=5 

so that, 

0sI sL-—1 
and 

Das 

2.B=L-1 

It can be proved that 
| j=1 
N = Ny 


almost surely. In pratice, we verified that Ny) is an 
increasing sequence for low N (N=10 for L= 500). 
Interpreting the increments 1,)J as "frequencies” in a 
multinomial experiment, we introduce the statistic: 


L, 


q = > 


where 


Ú 
H 


E=E 


This measure has the advantage of avoiding 
variance computations. 
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